#This program solves the Navier Stokes equation. It uses the method of
#constructed solutions to test for convergence, which should be 2nd order 

#from __future__ import division
import numpy as np
import numpy.fft as fft
from numpy import pi 
from pyfft.cl import Plan
import pyopencl as cl
import pyopencl.array as cl_array
from pyopencl.elementwise import ElementwiseKernel
import time
import sys
from pyvisfile.vtk import write_structured_grid
from pyvisfile.vtk import (
           UnstructuredGrid, DataArray,
           AppendedDataXMLGenerator,
           VTK_VERTEX, VF_LIST_OF_VECTORS, VF_LIST_OF_COMPONENTS)
from warnings import filterwarnings
filterwarnings("ignore", module="pyopencl.cache", lineno=336)
filterwarnings("ignore", module="pyopencl", lineno=163)
np.set_printoptions(threshold=sys.maxint)
TIMING = 0 
WGPSIZE_1 = 4 
WGPSIZE_2 = 4 
WGPSIZE_3 = 4 
WGPSIZE_11 = 4 
WGPSIZE_21 = 4
b_WGPSIZE = 2 
PR_INT = 80 
start_size = 5
end_size = 5 

def spread_routine(q, firstn1, nextn1, firstn2, nextn2, F_1, F_2, F_3, X_1, X_2, X_3, N, dx, dtheta, f_1, f_2, f_3):
    for shift0 in range(0,4):
        for shift1 in range(0,4):
            spread_kernel(q, (N/4, N/4), (WGPSIZE_11, WGPSIZE_21), firstn1.data, nextn1.data, firstn2.data, nextn2.data, F_1.data, F_2.data, F_3.data, X_1.data, X_2.data, X_3.data, N, dx, dtheta, shift0, shift1, f_1.data, f_2.data, f_3.data)

def update_routine(q, firstn1, nextn1, firstn2, nextn2, X_1, X_2, N, dx):
    for shift0 in range(0,4):
        for shift1 in range(0,4):
            updaten_kernel(q, (N/4, N/4), (WGPSIZE_11, WGPSIZE_21), firstn1.data, nextn1.data, firstn2.data, nextn2.data, X_1.data, X_2.data, N, dx, shift0, shift1)

def calc_transf(dx):
    N = np.int(np.round(1/dx,0))
    k = np.arange(0, N/2+1)
    k = np.hstack([k, np.arange(-N/2+1, 0)])
    k1 = np.tile(k, (N, N, 1))
    k2 = np.tile(k.reshape(N,1), (N, 1, N))
    k3 = np.tile(k.reshape(N,1,1), (1, N, N))
    L_tf = -4./dx**2 * (np.sin(pi*k1/N)**2+np.sin(pi*k2/N)**2+np.sin(pi*k3/N)**2) 
    L_tf[0,0,0] = 1. #avoid dividing by zero
    Nyq_filt = np.ones((N,N,N))
    Nyq_filt[N/2,:,:]=0
    Nyq_filt[:,N/2,:]=0 
    Nyq_filt[:,:,N/2]=0 
    p_filt = np.ones((N,N,N))
    p_filt[0,0,0] = 0
    return L_tf, Nyq_filt, p_filt

def solve(t_step):
    global dx, dt, mu, rho, b_q, dtheta
    N = np.int(np.round(1/dx,0))
 
    imag_gpu(u_1_gpu, u_1_tf_gpu)
    imag_gpu(u_2_gpu, u_2_tf_gpu)
    imag_gpu(u_3_gpu, u_3_tf_gpu)

    plan.execute(u_1_tf_gpu.data)
    plan.execute(u_2_tf_gpu.data)
    plan.execute(u_3_tf_gpu.data)
    
    calc_lagF_gpu(X_1_gpu, X_2_gpu, rr_lag_gpu, ip_s_lag_gpu, im_s_lag_gpu, ds, F_1_gpu, F_2_gpu)
    zero_gpu(F_3_gpu)

    interp_kernel(queue, (b_q, 1), (b_WGPSIZE, 1), u_1_gpu.data, u_2_gpu.data, u_3_gpu.data, X_1_gpu.data, X_2_gpu.data, X_3_gpu.data, N, dx, U_1_gpu.data, U_2_gpu.data, U_3_gpu.data)

    add_gpu(X_1_gpu, U_1_gpu, dt, X_1_tilde_gpu)
    add_gpu(X_2_gpu, U_2_gpu, dt, X_2_tilde_gpu)
    add_gpu(X_3_gpu, U_3_gpu, dt, X_3_tilde_gpu)

    copy_gpu(u_1_gpu,test_1_gpu) #for debugging 
    copy_gpu(u_2_gpu,test_2_gpu) 
    copy_gpu(u_3_gpu,test_3_gpu) 

    ss_conv_kernel(queue, (N,N,N), (WGPSIZE_1, WGPSIZE_2, WGPSIZE_3), u_1_gpu.data, u_2_gpu.data, u_3_gpu.data, dx, N, u_1_tilde_gpu.data, u_2_tilde_gpu.data, u_3_tilde_gpu.data)
    zero_gpu(f_1_gpu)
    zero_gpu(f_2_gpu)
    zero_gpu(f_3_gpu)
    zero_gpu(u_1_gpu)
    zero_gpu(u_2_gpu)
    zero_gpu(u_3_gpu)

    spread_routine(queue, firstn1_gpu, nextn1_gpu, firstn2_gpu, nextn2_gpu, F_1_gpu, F_2_gpu, F_3_gpu, X_1_gpu, X_2_gpu, X_3_gpu, N, dx, dtheta, u_1_gpu, u_2_gpu, u_3_gpu)
    #check power identity
    print 'power is', (U_1_gpu.get()*F_1_gpu.get()+U_2_gpu.get()*F_2_gpu.get()+U_3_gpu.get()*F_3_gpu.get()).sum()*dtheta-(test_1_gpu.get()*u_1_gpu.get()+test_2_gpu.get()*u_2_gpu.get()+test_3_gpu.get()*u_3_gpu.get()).sum()*dx*dx*dx
    #check momentum conservation
    print 'mom1 is', test_1_gpu.get().sum()*dx*dx*dx
    print 'mom2 is', test_2_gpu.get().sum()*dx*dx*dx
    print 'mom3 is', test_3_gpu.get().sum()*dx*dx*dx

    add_gpu(f_1_gpu, u_1_gpu, 1, f_1_gpu) 
    add_gpu(f_2_gpu, u_2_gpu, 1, f_2_gpu) 
    add_gpu(f_3_gpu, u_3_gpu, 1, f_3_gpu) 
    ##########################################################solve for u_tilde
    add_gpu(f_1_gpu, u_1_tilde_gpu, -rho, f_1_gpu)
    add_gpu(f_2_gpu, u_2_tilde_gpu, -rho, f_2_gpu)
    add_gpu(f_3_gpu, u_3_tilde_gpu, -rho, f_3_gpu)
    
    add_gpu(u_1_gpu, u_1_tilde_gpu, -rho, u_1_tilde_gpu)
    add_gpu(u_2_gpu, u_2_tilde_gpu, -rho, u_2_tilde_gpu)
    add_gpu(u_3_gpu, u_3_tilde_gpu, -rho, u_3_tilde_gpu)

    div_kernel(queue, (N,N,N), (WGPSIZE_1, WGPSIZE_2, WGPSIZE_3), f_1_gpu.data, f_2_gpu.data, f_3_gpu.data, dx, N, p_gpu.data)
    imag_gpu(p_gpu, p_tf_gpu)
    plan.execute(p_tf_gpu.data)

    L_adj_2_kernel(queue, (1,1), (1,1), L_tf_gpu.data)
    calc_p_gpu(p_tf_gpu, p_filt_gpu, L_tf_gpu, p_tf_gpu)
    plan.execute(p_tf_gpu.data, inverse=True)
    real_gpu(p_tf_gpu, p_gpu)

    grad_kernel(queue, (N,N,N), (WGPSIZE_1, WGPSIZE_2, WGPSIZE_3), p_gpu.data, dx, N, u_1_gpu.data, u_2_gpu.data, u_3_gpu.data)

    add_gpu(f_1_gpu, u_1_gpu, -1, f_1_gpu)
    add_gpu(f_2_gpu, u_2_gpu, -1, f_2_gpu)
    add_gpu(f_3_gpu, u_3_gpu, -1, f_3_gpu)

    imag_gpu(f_1_gpu, f_1_tf_gpu)
    imag_gpu(f_2_gpu, f_2_tf_gpu)
    imag_gpu(f_3_gpu, f_3_tf_gpu)

    plan.execute(f_1_tf_gpu.data)
    plan.execute(f_2_tf_gpu.data)
    plan.execute(f_3_tf_gpu.data)

    L_adj_1_kernel(queue, (1,1), (1,1), L_tf_gpu.data)
    solve_kernel(queue, (N,N,N), (WGPSIZE_1, WGPSIZE_2, WGPSIZE_3), u_1_tf_gpu.data, u_2_tf_gpu.data, u_3_tf_gpu.data, L_tf_gpu.data, f_1_tf_gpu.data, f_2_tf_gpu.data, f_3_tf_gpu.data, dt, rho, mu, N, u_1_tf_tilde_gpu.data, u_2_tf_tilde_gpu.data, u_3_tf_tilde_gpu.data)

    Nyq_gpu(u_1_tf_tilde_gpu, Nyq_filt_gpu, u_1_tf_tilde_gpu)
    Nyq_gpu(u_2_tf_tilde_gpu, Nyq_filt_gpu, u_2_tf_tilde_gpu)
    Nyq_gpu(u_3_tf_tilde_gpu, Nyq_filt_gpu, u_3_tf_tilde_gpu)

    zero_gpu(f_1_gpu)
    zero_gpu(f_2_gpu)
    zero_gpu(f_3_gpu)

    add_gpu(f_1_gpu, u_1_tilde_gpu, 0.5, f_1_gpu)
    add_gpu(f_2_gpu, u_2_tilde_gpu, 0.5, f_2_gpu)
    add_gpu(f_3_gpu, u_3_tilde_gpu, 0.5, f_3_gpu)

    plan.execute(u_1_tf_tilde_gpu.data, inverse=True)
    plan.execute(u_2_tf_tilde_gpu.data, inverse=True)
    plan.execute(u_3_tf_tilde_gpu.data, inverse=True)
    
    real_gpu(u_1_tf_tilde_gpu, u_1_tilde_gpu)
    real_gpu(u_2_tf_tilde_gpu, u_2_tilde_gpu)
    real_gpu(u_3_tf_tilde_gpu, u_3_tilde_gpu)
    ############################################################

    ############################################################solve for u_soln
    calc_lagF_gpu(X_1_tilde_gpu, X_2_tilde_gpu, rr_lag_gpu, ip_s_lag_gpu, im_s_lag_gpu, ds, F_1_gpu, F_2_gpu)
    zero_gpu(F_3_gpu)

    update_routine(queue, firstn1_gpu, nextn1_gpu, firstn2_gpu, nextn2_gpu, X_1_tilde_gpu, X_2_tilde_gpu, N, dx)

    zero_gpu(u_1_gpu)
    zero_gpu(u_2_gpu)
    zero_gpu(u_3_gpu)

    spread_routine(queue, firstn1_gpu, nextn1_gpu, firstn2_gpu, nextn2_gpu, F_1_gpu, F_2_gpu, F_3_gpu, X_1_tilde_gpu, X_2_tilde_gpu, X_3_tilde_gpu, N, dx, dtheta, u_1_gpu, u_2_gpu, u_3_gpu)

    add_gpu(f_1_gpu, u_1_gpu, 0.5, f_1_gpu)
    add_gpu(f_2_gpu, u_2_gpu, 0.5, f_2_gpu)
    add_gpu(f_3_gpu, u_3_gpu, 0.5, f_3_gpu)

    ss_conv_kernel(queue, (N,N,N), (WGPSIZE_1, WGPSIZE_2, WGPSIZE_3), u_1_tilde_gpu.data, u_2_tilde_gpu.data, u_3_tilde_gpu.data, dx, N, u_1_gpu.data, u_2_gpu.data, u_3_gpu.data)

    add_gpu(f_1_gpu, u_1_gpu, -rho*0.5, f_1_gpu)
    add_gpu(f_2_gpu, u_2_gpu, -rho*0.5, f_2_gpu)
    add_gpu(f_3_gpu, u_3_gpu, -rho*0.5, f_3_gpu)

    div_kernel(queue, (N,N,N), (WGPSIZE_1, WGPSIZE_2, WGPSIZE_3), f_1_gpu.data, f_2_gpu.data, f_3_gpu.data, dx, N, p_gpu.data)
    imag_gpu(p_gpu, p_tf_gpu)
    plan.execute(p_tf_gpu.data)

    L_adj_2_kernel(queue, (1,1), (1,1), L_tf_gpu.data)
    calc_p_gpu(p_tf_gpu, p_filt_gpu, L_tf_gpu, p_tf_gpu)
    plan.execute(p_tf_gpu.data, inverse=True)
    real_gpu(p_tf_gpu, p_gpu)

    grad_kernel(queue, (N,N,N), (WGPSIZE_1, WGPSIZE_2, WGPSIZE_3), p_gpu.data, dx, N, u_1_gpu.data, u_2_gpu.data, u_3_gpu.data)

    add_gpu(f_1_gpu, u_1_gpu, -1, f_1_gpu)
    add_gpu(f_2_gpu, u_2_gpu, -1, f_2_gpu)
    add_gpu(f_3_gpu, u_3_gpu, -1, f_3_gpu)

    imag_gpu(f_1_gpu, f_1_tf_gpu)
    imag_gpu(f_2_gpu, f_2_tf_gpu)
    imag_gpu(f_3_gpu, f_3_tf_gpu)

    plan.execute(f_1_tf_gpu.data)
    plan.execute(f_2_tf_gpu.data)
    plan.execute(f_3_tf_gpu.data)

    L_adj_1_kernel(queue, (1,1), (1,1), L_tf_gpu.data)
    solve_kernel(queue, (N,N,N), (WGPSIZE_1, WGPSIZE_2, WGPSIZE_3), u_1_tf_gpu.data, u_2_tf_gpu.data, u_3_tf_gpu.data, L_tf_gpu.data, f_1_tf_gpu.data, f_2_tf_gpu.data, f_3_tf_gpu.data, dt, rho, mu, N, u_1_tf_tilde_gpu.data, u_2_tf_tilde_gpu.data, u_3_tf_tilde_gpu.data)

    Nyq_gpu(u_1_tf_tilde_gpu, Nyq_filt_gpu, u_1_tf_tilde_gpu)
    Nyq_gpu(u_2_tf_tilde_gpu, Nyq_filt_gpu, u_2_tf_tilde_gpu)
    Nyq_gpu(u_3_tf_tilde_gpu, Nyq_filt_gpu, u_3_tf_tilde_gpu)

    plan.execute(u_1_tf_tilde_gpu.data, inverse=True)
    plan.execute(u_2_tf_tilde_gpu.data, inverse=True)
    plan.execute(u_3_tf_tilde_gpu.data, inverse=True)
    
    real_gpu(u_1_tf_tilde_gpu, u_1_gpu)
    real_gpu(u_2_tf_tilde_gpu, u_2_gpu)
    real_gpu(u_3_tf_tilde_gpu, u_3_gpu)

    interp_kernel(queue, (b_q, 1), (b_WGPSIZE, 1), u_1_gpu.data, u_2_gpu.data, u_3_gpu.data, X_1_tilde_gpu.data, X_2_tilde_gpu.data, X_3_tilde_gpu.data, N, dx, U_1_soln_gpu.data, U_2_soln_gpu.data, U_3_soln_gpu.data)
    add_gpu(U_1_gpu, U_1_soln_gpu, 1, U_1_gpu)
    add_gpu(U_2_gpu, U_2_soln_gpu, 1, U_2_gpu)
    add_gpu(U_3_gpu, U_3_soln_gpu, 1, U_3_gpu)
    add_gpu(X_1_gpu, U_1_gpu, dt/2., X_1_gpu)
    add_gpu(X_2_gpu, U_2_gpu, dt/2., X_2_gpu)
    add_gpu(X_3_gpu, U_3_gpu, dt/2., X_3_gpu)

    update_routine(queue, firstn1_gpu, nextn1_gpu, firstn2_gpu, nextn2_gpu, X_1_gpu, X_2_gpu, N, dx)
    ##############################################################
    return

def err(u_1, u_2, u_3, u_1_o, u_2_o, u_3_o):
    global dx
    NN_o = u_1_o.shape[0]
    errinf1 = 0.
    errinf2 = 0.
    errinf3 = 0.
    err11 = 0.
    err12 = 0.
    err13 = 0.
    err21 = 0.
    err22 = 0.
    err23 = 0.

    for i in range(0,NN_o):
        for j in range(0,NN_o):
            for k in range(0,NN_o):
                if (np.abs(u_1_o[k,i,j]-0.25*(u_1[2*k,2*i,2*j]+u_1[2*k,2*i+1,2*j]+u_1[2*k+1,2*i,2*j]+u_1[2*k+1,2*i+1,2*j]))>errinf1):
                    errinf1 = np.abs(u_1_o[k,i,j]-0.25*(u_1[2*k,2*i,2*j]+u_1[2*k,2*i+1,2*j]+u_1[2*k+1,2*i,2*j]+u_1[2*k+1,2*i+1,2*j])) 
                if (np.abs(u_2_o[k,i,j]-0.25*(u_2[2*k,2*i,2*j]+u_2[2*k,2*i,2*j+1]+u_2[2*k+1,2*i,2*j]+u_2[2*k+1,2*i,2*j+1]))>errinf2):
                    errinf2 = np.abs(u_2_o[k,i,j]-0.25*(u_2[2*k,2*i,2*j]+u_2[2*k,2*i,2*j+1]+u_2[2*k+1,2*i,2*j]+u_2[2*k+1,2*i,2*j+1])) 
                if (np.abs(u_3_o[k,i,j]-0.25*(u_3[2*k,2*i,2*j]+u_3[2*k,2*i+1,2*j]+u_3[2*k,2*i,2*j+1]+u_3[2*k,2*i+1,2*j+1]))>errinf3):
                    errinf3 = np.abs(u_3_o[k,i,j]-0.25*(u_3[2*k,2*i,2*j]+u_3[2*k,2*i+1,2*j]+u_3[2*k,2*i,2*j+1]+u_3[2*k,2*i+1,2*j+1])) 
                err21 += np.abs(u_1_o[k,i,j]-0.25*(u_1[2*k,2*i,2*j]+u_1[2*k,2*i+1,2*j]+u_1[2*k+1,2*i,2*j]+u_1[2*k+1,2*i+1,2*j]))**2 
                err22 += np.abs(u_2_o[k,i,j]-0.25*(u_2[2*k,2*i,2*j]+u_2[2*k,2*i,2*j+1]+u_2[2*k+1,2*i,2*j]+u_2[2*k+1,2*i,2*j+1]))**2 
                err23 += np.abs(u_3_o[k,i,j]-0.25*(u_3[2*k,2*i,2*j]+u_3[2*k,2*i+1,2*j]+u_3[2*k,2*i,2*j+1]+u_3[2*k,2*i+1,2*j+1]))**2 
                err11 += np.abs(u_1_o[k,i,j]-0.25*(u_1[2*k,2*i,2*j]+u_1[2*k,2*i+1,2*j]+u_1[2*k+1,2*i,2*j]+u_1[2*k+1,2*i+1,2*j])) 
                err12 += np.abs(u_2_o[k,i,j]-0.25*(u_2[2*k,2*i,2*j]+u_2[2*k,2*i,2*j+1]+u_2[2*k+1,2*i,2*j]+u_2[2*k+1,2*i,2*j+1])) 
                err13 += np.abs(u_3_o[k,i,j]-0.25*(u_3[2*k,2*i,2*j]+u_3[2*k,2*i+1,2*j]+u_3[2*k,2*i,2*j+1]+u_3[2*k,2*i+1,2*j+1])) 

    err21 = np.sqrt(err21*dx**3/8)
    err22 = np.sqrt(err22*dx**3/8)
    err23 = np.sqrt(err23*dx**3/8)
    err11 = (dx**3/8)*err11
    err12 = (dx**3/8)*err12
    err13 = (dx**3/8)*err13

    return errinf1, errinf2, errinf3, err11, err12, err13, err21, err22, err23 

def err_ex(clock_max, u_1, u_2, u_3):
    global dx, dt
    x_stag = np.arange(-0.5, 0.5, dx)
    xx_1_stag = np.tile(x_stag, (N, N, 1))
    xx_2_stag = np.tile(x_stag.reshape(N,1), (N, 1, N))
    xx_3_stag = np.tile(x_stag.reshape(N,1,1), (1, N, N))

    x_c = np.arange(-0.5+dx/2, 0.5+dx/2, dx)
    xx_1_c = np.tile(x_c, (N, N, 1))
    xx_2_c = np.tile(x_c.reshape(N,1), (N, 1, N))
    xx_3_c = np.tile(x_c.reshape(N,1,1), (1, N, N))
    
    u_1_ex = 2*np.cos(2*pi*xx_1_stag)*np.sin(4*pi*xx_2_c)*np.sin(clock_max*dt)*np.sin(2*pi*xx_3_c)
    u_2_ex = -np.sin(2*pi*xx_1_c)*np.cos(4*pi*xx_2_stag)*np.sin(clock_max*dt)*np.sin(2*pi*xx_3_c)
    u_3_ex = np.cos(2*pi*xx_1_c)*np.sin(4*pi*xx_2_c)*np.cos(clock_max*dt)

    errinf1 = np.abs(u_1-u_1_ex).max()
    errinf2 = np.abs(u_2-u_2_ex).max()
    errinf3 = np.abs(u_3-u_3_ex).max()
    err11 = np.sum(np.abs(u_1-u_1_ex))*dx**3
    err12 = np.sum(np.abs(u_2-u_2_ex))*dx**3
    err13 = np.sum(np.abs(u_3-u_3_ex))*dx**3
    err21 = np.sqrt(np.sum(np.abs(u_1-u_1_ex)**2*dx**3))
    err22 = np.sqrt(np.sum(np.abs(u_2-u_2_ex)**2*dx**3))
    err23 = np.sqrt(np.sum(np.abs(u_3-u_3_ex)**2*dx**3))

    return errinf1, errinf2, errinf3, err11, err12, err13, err21, err22, err23

ctx = cl.create_some_context()
queue = cl.CommandQueue(ctx)
Nyq_gpu = ElementwiseKernel(ctx,
	"float2 *u, float *filt, float2 *f",
	"f[i] = filt[i]*u[i]",
	"Nyq_gpu")

calc_lagF_gpu = ElementwiseKernel(ctx,
	"float *X_1, float *X_2, "
	"float *rr_lag, unsigned *ip_s_lag, unsigned *im_s_lag, " 
	"float ds, "
	"float *F_1, float *F_2",
	"F_1[i] = 1./(ds*ds)*(1.+sin(2.*pi*rr_lag[i]-pi/2.))"
	   "*(X_1[ip_s_lag[i]]-2.*X_1[i]+X_1[im_s_lag[i]]); "
	"F_2[i] = 1./(ds*ds)*(1.+sin(2.*pi*rr_lag[i]-pi/2.))"
	   "*(X_2[ip_s_lag[i]]-2.*X_2[i]+X_2[im_s_lag[i]])",
	#"F_1[i] = 1./(ds*ds)*(X_1[ip_s_lag[i]]-2.*X_1[i]+X_1[im_s_lag[i]]); "
	#"F_2[i] = 1./(ds*ds)*(X_2[ip_s_lag[i]]-2.*X_2[i]+X_2[im_s_lag[i]]) ",
        #"F_1[i] = 0.;"
        #"F_2[i] = 0.",
	"calc_lagF_gpu",
	preamble="""
	#define pi (float)3.141592653589793
	""")

calc_f_gpu = ElementwiseKernel(ctx,
        "float t, float dx, float mu, float rho, "
        "float *xx_1_stag, float *xx_2_stag, float *xx_3_stag,"
        "float *xx_1_c, float *xx_2_c, float *xx_3_c,"
        "float *f_1, float *f_2, float *f_3",
        "f_1[i] = rho*2*cos(2*pi*xx_1_stag[i])*sin(4*pi*xx_2_c[i])*sin(2*pi*xx_3_c[i])*cos(t)"
           "+2*pi*cos(2*pi*xx_1_stag[i])*sin(2*pi*xx_2_c[i])*cos(2*pi*xx_3_c[i])*cos(t)"
           "+mu*48*pi*pi*cos(2*pi*xx_1_stag[i])*sin(4*pi*xx_2_c[i])*sin(2*pi*xx_3_c[i])*sin(t)"
           "-rho*4*pi*sin(t)*sin(t)*sin(4*pi*xx_1_stag[i])*sin(2*pi*xx_3_c[i])*sin(2*pi*xx_3_c[i])"
           "+rho*4*pi*sin(t)*cos(t)*cos(2*pi*xx_1_stag[i])*cos(2*pi*xx_1_stag[i])*sin(4*pi*xx_2_c[i])"
           "*sin(4*pi*xx_2_c[i])*cos(2*pi*xx_3_c[i]);"
        "f_2[i] = -rho*sin(2*pi*xx_1_c[i])*cos(4*pi*xx_2_stag[i])*sin(2*pi*xx_3_c[i])*cos(t)"
           "+2*pi*sin(2*pi*xx_1_c[i])*cos(2*pi*xx_2_stag[i])*cos(2*pi*xx_3_c[i])*cos(t)"
           "-mu*24*pi*pi*sin(2*pi*xx_1_c[i])*cos(4*pi*xx_2_stag[i])*sin(2*pi*xx_3_c[i])*sin(t)"
           "-rho*2*pi*sin(t)*sin(t)*sin(8*pi*xx_2_stag[i])*sin(2*pi*xx_3_c[i])*sin(2*pi*xx_3_c[i])"
           "-rho*0.5*pi*sin(t)*cos(t)*sin(4*pi*xx_1_c[i])*sin(8*pi*xx_2_stag[i])"
           "*cos(2*pi*xx_3_c[i]);"
        "f_3[i] = -rho*cos(2*pi*xx_1_c[i])*sin(4*pi*xx_2_c[i])*sin(t)"
           "-2*pi*sin(2*pi*xx_1_c[i])*sin(2*pi*xx_2_c[i])*sin(2*pi*xx_3_stag[i])*cos(t)"
           "+mu*20*pi*pi*cos(2*pi*xx_1_c[i])*sin(4*pi*xx_2_c[i])*cos(t)"
           "-rho*2*pi*cos(t)*sin(t)*sin(4*pi*xx_1_c[i])*sin(2*pi*xx_3_stag[i])",
        "calc_f_gpu",
        preamble="""
        #define pi (float)3.141592653589793
        """)

complex_mult_gpu = ElementwiseKernel(ctx,
	"float2 *u_1, float2 *u_2, float2 *f",
	"f[i] = complex_mul(u_1[i], u_2[i])",
	"complex_mult_gpu",
	preamble="""
	#define complex_ctr(x, y) (float2)(x, y)
	#define complex_mul(a, b) complex_ctr(-(a).y*(b).y+(a).x*(b).x, (a).y*(b).x+(a).x*(b).y)
	""")

mult_gpu = ElementwiseKernel(ctx,
	"float *u_1, float *u_2, float *f",
	"f[i] = u_1[i]*u_2[i]",
	"mult_gpu")

add_gpu = ElementwiseKernel(ctx,
	"float *u_1, float *u_2, float a, float *f",
	"f[i] = u_1[i]+a*u_2[i]",
	"add_gpu")

copy_gpu = ElementwiseKernel(ctx,
	"float *u, float *f",
	"f[i] = u[i]",
	"copy_gpu")
    
real_gpu = ElementwiseKernel(ctx,
	"float2 *x, float *z",
	"z[i] = x[i].x",
	"real_gpu")

imag_gpu = ElementwiseKernel(ctx,
	"float *x, float2 *z",
	"z[i].x = x[i];"
	"z[i].y = 0",
	"imag_gpu")

zero_gpu = ElementwiseKernel(ctx,
	"float *x",
	"x[i] = 0",
	"zero_gpu")
 
null_imag_gpu = ElementwiseKernel(ctx,
	"float2 *x, float2 *z",
	"z[i].x = x[i].x;"
	"z[i].y = 0",
	"null_imag_gpu")

calc_p_gpu = ElementwiseKernel(ctx,
	"float2 *w_tf, float *filt, float2 *L_tf, float2 *p_tf",
	"p_tf[i] = filt[i]*complex_div(w_tf[i], L_tf[i])",
	"calc_p_gpu",
	preamble="""
	#define complex_ctr(x, y) (float2)(x, y)
	#define complex_div(a, b) 1/((b).x*(b).x+(b).y*(b).y)*complex_ctr((a).y*(b).y+(a).x*(b).x,(a).y*(b).x-(a).x*(b).y)
	""")

solve_prg = cl.Program(ctx, """
    //#pragma OPENCL EXTENSION cl_amd_printf : enable
    #define complex_ctr(x, y) (float2)(x, y)
    #define complex_mul(a, b) complex_ctr(-(a).y*(b).y+(a).x*(b).x, (a).y*(b).x+(a).x*(b).y)
    #define complex_div(a, b) 1/((b).x*(b).x+(b).y*(b).y)*complex_ctr((a).y*(b).y+(a).x*(b).x,(a).y*(b).x-(a).x*(b).y)
    __kernel void solve(__global const float2 *u_1_tf, __global const float2 *u_2_tf, 
    __global const float2 *u_3_tf, __global const float2 *L_tf,
    __global const float2 *f_1_tf, __global const float2 *f_2_tf, 
    __global const float2 *f_3_tf, const float dt, 
    const float rho, const float mu, const unsigned N,
    __global float2 *u_1_tf_tilde, __global float2 *u_2_tf_tilde, __global float2 *u_3_tf_tilde)
    {
      unsigned i = get_global_id(0);
      unsigned j = get_global_id(1);
      unsigned k = get_global_id(2);
      unsigned base = i + N*j + N*N*k;
      float2 m;
      m.x = 1.0+mu*dt/(2.0*rho)*L_tf[base].x; 
      m.y = mu*dt/(2.0*rho)*L_tf[base].y; 
      float2 rhs_1 = complex_mul(m, u_1_tf[base])+dt/rho*f_1_tf[base];
      float2 rhs_2 = complex_mul(m, u_2_tf[base])+dt/rho*f_2_tf[base];
      float2 rhs_3 = complex_mul(m, u_3_tf[base])+dt/rho*f_3_tf[base];
      float2 lhs;
      lhs.x = 1.0-mu*dt/(2.0*rho)*L_tf[base].x;
      lhs.y = -mu*dt/(2.0*rho)*L_tf[base].y;
      u_1_tf_tilde[base] = complex_div(rhs_1, lhs);
      u_2_tf_tilde[base] = complex_div(rhs_2, lhs);
      u_3_tf_tilde[base] = complex_div(rhs_3, lhs);
    }
    """).build()

prg = cl.Program(ctx, """
    //#pragma OPENCL EXTENSION cl_amd_printf : enable
    __kernel void div(__global const float *u_1,
    __global const float *u_2, __global const float *u_3, 
    const float dx, const unsigned N, __global float *f)
    {
      unsigned i = get_global_id(0);
      unsigned j = get_global_id(1);
      unsigned k = get_global_id(2);
      //f[i,j] = 1/dx*(u_1[i_2[i,j],ip_1[i,j]]-u_1[i,j]+u_2[ip_2[i,j],i_1[i,j]]-u_2[i,j];
      unsigned base = i + N*j + N*N*k;
      unsigned base1 = (i+1 & N-1)+N*j + N*N*k;
      unsigned base2 = i+N*(j+1 & N-1) + N*N*k;
      unsigned base3 = i + N*j + N*N*(k+1 & N-1);
      f[base] = 1/dx*(u_1[base1]-u_1[base]+u_2[base2]-u_2[base]+u_3[base3]-u_3[base]);
    }
    __kernel void L_adj_1(__global float2 *L_tf)
    {
      L_tf[0].x = 0.0;
      L_tf[0].y = 0.0;
    }
    __kernel void L_adj_2(__global float2 *L_tf)
    {
      L_tf[0].x = 1.0;
      L_tf[0].y = 0.0;
    }
    __kernel void grad(__global const float *phi,
    const float dx, const unsigned N,
    __global float *g_1, __global float *g_2, __global float *g_3)
    {
      unsigned i = get_global_id(0);
      unsigned j = get_global_id(1);
      unsigned k = get_global_id(2);
      //grad_1 = 1/dx*(phi[i_2,i_1]-phi[i_2,im_1])
      //grad_2 = 1/dx*(phi[i_2,i_1]-phi[im_2,i_1])
      unsigned base = i + N*j + N*N*k;
      unsigned base1 = (i-1 & N-1)+N*j + N*N*k;
      unsigned base2 = i+N*(j-1 & N-1) + N*N*k;
      unsigned base3 = i + N*j + N*N*(k-1 & N-1);
      g_1[base] = 1/dx*(phi[base]-phi[base1]);
      g_2[base] = 1/dx*(phi[base]-phi[base2]);
      g_3[base] = 1/dx*(phi[base]-phi[base3]);
    }
    __kernel void ss_conv(__global const float *u_1,
    __global const float *u_2, __global const float *u_3,
    const float dx, const unsigned N,
    __global float *S_1, __global float *S_2, __global float *S_3)
    {
      unsigned i = get_global_id(0);
      unsigned j = get_global_id(1);
      unsigned k = get_global_id(2);
      //Compute S = 0.5*(div(UU) + (U*grad)U).
      unsigned base = i+N*j+N*N*k;
      unsigned base1 = (i-1 & N-1)+N*j+N*N*k;
      unsigned base2 = (i+1 & N-1)+N*j+N*N*k;
      unsigned base3 = i+N*(j-1 & N-1)+N*N*k;
      unsigned base4 = i+N*(j+1 & N-1)+N*N*k;
      unsigned base5 = (i-1 & N-1)+N*(j+1 & N-1)+N*N*k;
      unsigned base6 = (i+1 & N-1)+N*(j-1 & N-1)+N*N*k;
      unsigned base7 = i+N*j+N*N*(k-1 & N-1);
      unsigned base8 = i+N*j+N*N*(k+1 & N-1);
      unsigned base9 = (i-1 & N-1)+N*j+N*N*(k+1 & N-1);
      unsigned base10 = i+N*(j-1 & N-1)+N*N*(k+1 & N-1);
      unsigned base11 = (i+1 & N-1)+N*j+N*N*(k-1 & N-1);
      unsigned base12 = i+N*(j+1 & N-1)+N*N*(k-1 & N-1);
      float D_x, D_y, D_z, A_x, A_y, A_z, D, A;
      D_x = (0.25/dx)*(
          (u_1[base2]+u_1[base])*(u_1[base2]+u_1[base]) -
          (u_1[base]+u_1[base1])*(u_1[base]+u_1[base1]));
      D_y = (0.25/dx)*(
          (u_2[base4]+u_2[base5])*(u_1[base4]+u_1[base]) -
          (u_2[base]+u_2[base1])*(u_1[base]+u_1[base3]));
      D_z = (0.25/dx)*(
          (u_3[base8]+u_3[base9])*(u_1[base8]+u_1[base]) -
          (u_3[base]+u_3[base1])*(u_1[base]+u_1[base7]));
      D = D_x + D_y + D_z;
      A_x =  (0.25/dx)*(
          (u_1[base2]+u_1[base])*(u_1[base2]-u_1[base]) +
          (u_1[base]+u_1[base1])*(u_1[base]-u_1[base1]));
      A_y = (0.25/dx)*(
            (u_2[base4]+u_2[base5])*(u_1[base4]-u_1[base]) +
            (u_2[base]+u_2[base1])*(u_1[base]-u_1[base3]));
      A_z = (0.25/dx)*(
            (u_3[base8]+u_3[base9])*(u_1[base8]-u_1[base]) +
            (u_3[base]+u_3[base1])*(u_1[base]-u_1[base7]));
      A = A_x + A_y + A_z;
      S_1[base] = 0.5*(D+A);
      D_x = (0.25/dx)*(
          (u_1[base2]+u_1[base6])*(u_2[base2]+u_2[base]) -
          (u_1[base]+u_1[base3])*(u_2[base]+u_2[base1]));
      D_y = (0.25/dx)*(
          (u_2[base4]+u_2[base])*(u_2[base4]+u_2[base]) -
          (u_2[base]+u_2[base3])*(u_2[base]+u_2[base3]));
      D_z = (0.25/dx)*(
          (u_3[base8]+u_3[base10])*(u_2[base8]+u_2[base]) -
          (u_3[base]+u_3[base3])*(u_2[base]+u_2[base7]));
      D = D_x + D_y + D_z;
      A_x =  (0.25/dx)*(
          (u_1[base2]+u_1[base6])*(u_2[base2]-u_2[base]) +
          (u_1[base]+u_1[base3])*(u_2[base]-u_2[base1]));
      A_y = (0.25/dx)*(
            (u_2[base4]+u_2[base])*(u_2[base4]-u_2[base]) +
            (u_2[base]+u_2[base3])*(u_2[base]-u_2[base3]));
      A_z = (0.25/dx)*(
            (u_3[base8]+u_3[base10])*(u_2[base8]-u_2[base]) +
            (u_3[base]+u_3[base3])*(u_2[base]-u_2[base7]));
      A = A_x + A_y + A_z;
      S_2[base] = 0.5*(D+A);
      D_x = (0.25/dx)*(
          (u_1[base2]+u_1[base11])*(u_3[base2]+u_3[base]) -
          (u_1[base]+u_1[base7])*(u_3[base]+u_3[base1]));
      D_y = (0.25/dx)*(
          (u_2[base4]+u_2[base12])*(u_3[base4]+u_3[base]) -
          (u_2[base]+u_2[base7])*(u_3[base]+u_3[base3]));
      D_z = (0.25/dx)*(
          (u_3[base8]+u_3[base])*(u_3[base8]+u_3[base]) -
          (u_3[base]+u_3[base7])*(u_3[base]+u_3[base7]));
      D = D_x + D_y + D_z;
      A_x =  (0.25/dx)*(
          (u_1[base2]+u_1[base11])*(u_3[base2]-u_3[base]) +
          (u_1[base]+u_1[base7])*(u_3[base]-u_3[base1]));
      A_y = (0.25/dx)*(
            (u_2[base4]+u_2[base12])*(u_3[base4]-u_3[base]) +
            (u_2[base]+u_2[base7])*(u_3[base]-u_3[base3]));
      A_z = (0.25/dx)*(
            (u_3[base8]+u_3[base])*(u_3[base8]-u_3[base]) +
            (u_3[base]+u_3[base7])*(u_3[base]-u_3[base7]));
      A = A_x + A_y + A_z;
      S_3[base] = 0.5*(D+A);
    }
    """).build()

ib_prg = cl.Program(ctx, """
    //#pragma OPENCL EXTENSION cl_amd_printf : enable
    __kernel void spread(__global const int *firstn1, __global const int *nextn1, __global const int *firstn2, __global const int *nextn2, __global const float *F_1, __global const float *F_2, __global const float *F_3, __global const float *X_1, __global const float *X_2, __global const float *X_3, const unsigned N, const float dx, const float dtheta, const unsigned shift0, const unsigned shift1, __global float *f_1, __global float *f_2, __global float *f_3)
    {
      unsigned i = 4*get_global_id(0)+shift0;
      unsigned j = 4*get_global_id(1)+shift1;

      int k1 = firstn1[i+N*j]; 
      int k2 = firstn2[i+N*j]; 
      int k3 = -1; //for now

      if (k1 == -1 && k2 == -1 && k3 == -1){
      return;
      }

      float c = dtheta/(dx*dx*dx);

      int ind_1x, ind_1y, ind_1z, ind_1w;
      int ind_2x, ind_2y, ind_2z, ind_2w;
      int ind_3x, ind_3y, ind_3z, ind_3w;

      float q1, q2, q3;

      float v1x, v1y, v1z, v1w;
      float v2x, v2y, v2z, v2w;
      float v3x, v3y, v3z, v3w;

      float ws01, ws11, ws21, ws31;
      float ws41, ws51, ws61, ws71;
      float ws81, ws91, wsa1, wsb1;
      float wsc1, wsd1, wse1, wsf1;

      float ws02, ws12, ws22, ws32;
      float ws42, ws52, ws62, ws72;
      float ws82, ws92, wsa2, wsb2;
      float wsc2, wsd2, wse2, wsf2;

      float ws03, ws13, ws23, ws33;
      float ws43, ws53, ws63, ws73;
      float ws83, ws93, wsa3, wsb3;
      float wsc3, wsd3, wse3, wsf3;

      float ws04, ws14, ws24, ws34;
      float ws44, ws54, ws64, ws74;
      float ws84, ws94, wsa4, wsb4;
      float wsc4, wsd4, wse4, wsf4;

      while(k1 != -1)
      {
      float s_1_stag = (X_1[k1]+0.5)/dx;
      float s_2_stag = (X_2[k1]+0.5)/dx;       
      float s_3_stag = (X_3[k1]+0.5)/dx;       
      int i_1_stag = convert_int(floor(s_1_stag));
      int i_2_stag = convert_int(floor(s_2_stag));
      int i_3_stag = convert_int(floor(s_3_stag));
      float r_1_stag = s_1_stag-i_1_stag;
      float r_2_stag = s_2_stag-i_2_stag;
      float r_3_stag = s_3_stag-i_3_stag;
      float s_1_c = (X_1[k1]+0.5-dx/2)/dx;
      float s_2_c = (X_2[k1]+0.5-dx/2)/dx;        
      float s_3_c = (X_3[k1]+0.5-dx/2)/dx;        
      int i_1_c = convert_int(floor(s_1_c));
      int i_2_c = convert_int(floor(s_2_c));
      int i_3_c = convert_int(floor(s_3_c));
      float r_1_c = s_1_c-i_1_c;
      float r_2_c = s_2_c-i_2_c;
      float r_3_c = s_3_c-i_3_c;
      
      ind_1x = (i_1_stag-1 & N-1);
      ind_1y = (i_1_stag & N-1);
      ind_1z = (i_1_stag+1 & N-1);
      ind_1w = (i_1_stag+2 & N-1);

      ind_2x = (i_2_c-1 & N-1);
      ind_2y = (i_2_c & N-1);
      ind_2z = (i_2_c+1 & N-1);
      ind_2w = (i_2_c+2 & N-1);

      ind_3x = (i_3_c-1 & N-1);
      ind_3y = (i_3_c & N-1);
      ind_3z = (i_3_c+1 & N-1);
      ind_3w = (i_3_c+2 & N-1);

      q1 = sqrt(1.+4.0*r_1_stag*(1-r_1_stag)); 
      q2 = sqrt(1.+4.0*r_2_c*(1-r_2_c));
      q3 = sqrt(1.+4.0*r_3_c*(1-r_3_c));

      v1x = 3.-2.*r_1_stag-q1; 
      v1y = 3.-2.*r_1_stag+q1; 
      v1z = 1.+2.*r_1_stag+q1; 
      v1w = 1.+2.*r_1_stag-q1; 
      v1x = v1x/8.;
      v1y = v1y/8.;
      v1z = v1z/8.;
      v1w = v1w/8.;

      v2x = 3.-2.*r_2_c-q2; 
      v2y = 3.-2.*r_2_c+q2; 
      v2z = 1.+2.*r_2_c+q2; 
      v2w = 1.+2.*r_2_c-q2; 
      v2x = v2x/8.;
      v2y = v2y/8.;
      v2z = v2z/8.;
      v2w = v2w/8.;

      v3x = 3.-2.*r_3_c-q2; 
      v3y = 3.-2.*r_3_c+q2; 
      v3z = 1.+2.*r_3_c+q2; 
      v3w = 1.+2.*r_3_c-q2; 
      v3x = v3x/8.;
      v3y = v3y/8.;
      v3z = v3z/8.;
      v3w = v3w/8.;
      
      ws01 = v1x*v2x*v3x;
      ws11 = v1y*v2x*v3x;
      ws21 = v1z*v2x*v3x;
      ws31 = v1w*v2x*v3x;
      ws41 = v1x*v2y*v3x;
      ws51 = v1y*v2y*v3x;
      ws61 = v1z*v2y*v3x;
      ws71 = v1w*v2y*v3x;
      ws81 = v1x*v2z*v3x;
      ws91 = v1y*v2z*v3x;
      wsa1 = v1z*v2z*v3x;
      wsb1 = v1w*v2z*v3x;
      wsc1 = v1x*v2w*v3x;
      wsd1 = v1y*v2w*v3x;
      wse1 = v1z*v2w*v3x;
      wsf1 = v1w*v2w*v3x;
      
      ws02 = v1x*v2x*v3y;
      ws12 = v1y*v2x*v3y;
      ws22 = v1z*v2x*v3y;
      ws32 = v1w*v2x*v3y;
      ws42 = v1x*v2y*v3y;
      ws52 = v1y*v2y*v3y;
      ws62 = v1z*v2y*v3y;
      ws72 = v1w*v2y*v3y;
      ws82 = v1x*v2z*v3y;
      ws92 = v1y*v2z*v3y;
      wsa2 = v1z*v2z*v3y;
      wsb2 = v1w*v2z*v3y;
      wsc2 = v1x*v2w*v3y;
      wsd2 = v1y*v2w*v3y;
      wse2 = v1z*v2w*v3y;
      wsf2 = v1w*v2w*v3y;

      ws03 = v1x*v2x*v3z;
      ws13 = v1y*v2x*v3z;
      ws23 = v1z*v2x*v3z;
      ws33 = v1w*v2x*v3z;
      ws43 = v1x*v2y*v3z;
      ws53 = v1y*v2y*v3z;
      ws63 = v1z*v2y*v3z;
      ws73 = v1w*v2y*v3z;
      ws83 = v1x*v2z*v3z;
      ws93 = v1y*v2z*v3z;
      wsa3 = v1z*v2z*v3z;
      wsb3 = v1w*v2z*v3z;
      wsc3 = v1x*v2w*v3z;
      wsd3 = v1y*v2w*v3z;
      wse3 = v1z*v2w*v3z;
      wsf3 = v1w*v2w*v3z;

      ws04 = v1x*v2x*v3w;
      ws14 = v1y*v2x*v3w;
      ws24 = v1z*v2x*v3w;
      ws34 = v1w*v2x*v3w;
      ws44 = v1x*v2y*v3w;
      ws54 = v1y*v2y*v3w;
      ws64 = v1z*v2y*v3w;
      ws74 = v1w*v2y*v3w;
      ws84 = v1x*v2z*v3w;
      ws94 = v1y*v2z*v3w;
      wsa4 = v1z*v2z*v3w;
      wsb4 = v1w*v2z*v3w;
      wsc4 = v1x*v2w*v3w;
      wsd4 = v1y*v2w*v3w;
      wse4 = v1z*v2w*v3w;
      wsf4 = v1w*v2w*v3w;

      float my_F_1 = F_1[k1];
      f_1[ind_1x+N*ind_2x+N*N*ind_3x] += c*my_F_1*ws01;
      f_1[ind_1y+N*ind_2x+N*N*ind_3x] += c*my_F_1*ws11;
      f_1[ind_1z+N*ind_2x+N*N*ind_3x] += c*my_F_1*ws21;
      f_1[ind_1w+N*ind_2x+N*N*ind_3x] += c*my_F_1*ws31;
      f_1[ind_1x+N*ind_2y+N*N*ind_3x] += c*my_F_1*ws41;
      f_1[ind_1y+N*ind_2y+N*N*ind_3x] += c*my_F_1*ws51;
      f_1[ind_1z+N*ind_2y+N*N*ind_3x] += c*my_F_1*ws61;
      f_1[ind_1w+N*ind_2y+N*N*ind_3x] += c*my_F_1*ws71;
      f_1[ind_1x+N*ind_2z+N*N*ind_3x] += c*my_F_1*ws81;
      f_1[ind_1y+N*ind_2z+N*N*ind_3x] += c*my_F_1*ws91;
      f_1[ind_1z+N*ind_2z+N*N*ind_3x] += c*my_F_1*wsa1;
      f_1[ind_1w+N*ind_2z+N*N*ind_3x] += c*my_F_1*wsb1;
      f_1[ind_1x+N*ind_2w+N*N*ind_3x] += c*my_F_1*wsc1;
      f_1[ind_1y+N*ind_2w+N*N*ind_3x] += c*my_F_1*wsd1;
      f_1[ind_1z+N*ind_2w+N*N*ind_3x] += c*my_F_1*wse1;
      f_1[ind_1w+N*ind_2w+N*N*ind_3x] += c*my_F_1*wsf1;

      f_1[ind_1x+N*ind_2x+N*N*ind_3y] += c*my_F_1*ws02;
      f_1[ind_1y+N*ind_2x+N*N*ind_3y] += c*my_F_1*ws12;
      f_1[ind_1z+N*ind_2x+N*N*ind_3y] += c*my_F_1*ws22;
      f_1[ind_1w+N*ind_2x+N*N*ind_3y] += c*my_F_1*ws32;
      f_1[ind_1x+N*ind_2y+N*N*ind_3y] += c*my_F_1*ws42;
      f_1[ind_1y+N*ind_2y+N*N*ind_3y] += c*my_F_1*ws52;
      f_1[ind_1z+N*ind_2y+N*N*ind_3y] += c*my_F_1*ws62;
      f_1[ind_1w+N*ind_2y+N*N*ind_3y] += c*my_F_1*ws72;
      f_1[ind_1x+N*ind_2z+N*N*ind_3y] += c*my_F_1*ws82;
      f_1[ind_1y+N*ind_2z+N*N*ind_3y] += c*my_F_1*ws92;
      f_1[ind_1z+N*ind_2z+N*N*ind_3y] += c*my_F_1*wsa2;
      f_1[ind_1w+N*ind_2z+N*N*ind_3y] += c*my_F_1*wsb2;
      f_1[ind_1x+N*ind_2w+N*N*ind_3y] += c*my_F_1*wsc2;
      f_1[ind_1y+N*ind_2w+N*N*ind_3y] += c*my_F_1*wsd2;
      f_1[ind_1z+N*ind_2w+N*N*ind_3y] += c*my_F_1*wse2;
      f_1[ind_1w+N*ind_2w+N*N*ind_3y] += c*my_F_1*wsf2;

      f_1[ind_1x+N*ind_2x+N*N*ind_3z] += c*my_F_1*ws03;
      f_1[ind_1y+N*ind_2x+N*N*ind_3z] += c*my_F_1*ws13;
      f_1[ind_1z+N*ind_2x+N*N*ind_3z] += c*my_F_1*ws23;
      f_1[ind_1w+N*ind_2x+N*N*ind_3z] += c*my_F_1*ws33;
      f_1[ind_1x+N*ind_2y+N*N*ind_3z] += c*my_F_1*ws43;
      f_1[ind_1y+N*ind_2y+N*N*ind_3z] += c*my_F_1*ws53;
      f_1[ind_1z+N*ind_2y+N*N*ind_3z] += c*my_F_1*ws63;
      f_1[ind_1w+N*ind_2y+N*N*ind_3z] += c*my_F_1*ws73;
      f_1[ind_1x+N*ind_2z+N*N*ind_3z] += c*my_F_1*ws83;
      f_1[ind_1y+N*ind_2z+N*N*ind_3z] += c*my_F_1*ws93;
      f_1[ind_1z+N*ind_2z+N*N*ind_3z] += c*my_F_1*wsa3;
      f_1[ind_1w+N*ind_2z+N*N*ind_3z] += c*my_F_1*wsb3;
      f_1[ind_1x+N*ind_2w+N*N*ind_3z] += c*my_F_1*wsc3;
      f_1[ind_1y+N*ind_2w+N*N*ind_3z] += c*my_F_1*wsd3;
      f_1[ind_1z+N*ind_2w+N*N*ind_3z] += c*my_F_1*wse3;
      f_1[ind_1w+N*ind_2w+N*N*ind_3z] += c*my_F_1*wsf3;

      f_1[ind_1x+N*ind_2x+N*N*ind_3w] += c*my_F_1*ws04;
      f_1[ind_1y+N*ind_2x+N*N*ind_3w] += c*my_F_1*ws14;
      f_1[ind_1z+N*ind_2x+N*N*ind_3w] += c*my_F_1*ws24;
      f_1[ind_1w+N*ind_2x+N*N*ind_3w] += c*my_F_1*ws34;
      f_1[ind_1x+N*ind_2y+N*N*ind_3w] += c*my_F_1*ws44;
      f_1[ind_1y+N*ind_2y+N*N*ind_3w] += c*my_F_1*ws54;
      f_1[ind_1z+N*ind_2y+N*N*ind_3w] += c*my_F_1*ws64;
      f_1[ind_1w+N*ind_2y+N*N*ind_3w] += c*my_F_1*ws74;
      f_1[ind_1x+N*ind_2z+N*N*ind_3w] += c*my_F_1*ws84;
      f_1[ind_1y+N*ind_2z+N*N*ind_3w] += c*my_F_1*ws94;
      f_1[ind_1z+N*ind_2z+N*N*ind_3w] += c*my_F_1*wsa4;
      f_1[ind_1w+N*ind_2z+N*N*ind_3w] += c*my_F_1*wsb4;
      f_1[ind_1x+N*ind_2w+N*N*ind_3w] += c*my_F_1*wsc4;
      f_1[ind_1y+N*ind_2w+N*N*ind_3w] += c*my_F_1*wsd4;
      f_1[ind_1z+N*ind_2w+N*N*ind_3w] += c*my_F_1*wse4;
      f_1[ind_1w+N*ind_2w+N*N*ind_3w] += c*my_F_1*wsf4;

      k1 = nextn1[k1];
      }

      while(k2 != -1)
      {
      float s_1_stag = (X_1[k2]+0.5)/dx;
      float s_2_stag = (X_2[k2]+0.5)/dx;       
      float s_3_stag = (X_3[k2]+0.5)/dx;       
      int i_1_stag = convert_int(floor(s_1_stag));
      int i_2_stag = convert_int(floor(s_2_stag));
      int i_3_stag = convert_int(floor(s_3_stag));
      float r_1_stag = s_1_stag-i_1_stag;
      float r_2_stag = s_2_stag-i_2_stag;
      float r_3_stag = s_3_stag-i_3_stag;
      float s_1_c = (X_1[k2]+0.5-dx/2)/dx;
      float s_2_c = (X_2[k2]+0.5-dx/2)/dx;        
      float s_3_c = (X_3[k2]+0.5-dx/2)/dx;        
      int i_1_c = convert_int(floor(s_1_c));
      int i_2_c = convert_int(floor(s_2_c));
      int i_3_c = convert_int(floor(s_3_c));
      float r_1_c = s_1_c-i_1_c;
      float r_2_c = s_2_c-i_2_c;
      float r_3_c = s_3_c-i_3_c;
      
      ind_1x = (i_1_c-1 & N-1);
      ind_1y = (i_1_c & N-1);
      ind_1z = (i_1_c+1 & N-1);
      ind_1w = (i_1_c+2 & N-1);

      ind_2x = (i_2_stag-1 & N-1);
      ind_2y = (i_2_stag & N-1);
      ind_2z = (i_2_stag+1 & N-1);
      ind_2w = (i_2_stag+2 & N-1);

      ind_3x = (i_3_c-1 & N-1);
      ind_3y = (i_3_c & N-1);
      ind_3z = (i_3_c+1 & N-1);
      ind_3w = (i_3_c+2 & N-1);

      q1 = sqrt(1.+4.0*r_1_c*(1-r_1_c)); 
      q2 = sqrt(1.+4.0*r_2_stag*(1-r_2_stag));
      q3 = sqrt(1.+4.0*r_3_c*(1-r_3_c));

      v1x = 3.-2.*r_1_c-q1; 
      v1y = 3.-2.*r_1_c+q1; 
      v1z = 1.+2.*r_1_c+q1; 
      v1w = 1.+2.*r_1_c-q1; 
      v1x = v1x/8.;
      v1y = v1y/8.;
      v1z = v1z/8.;
      v1w = v1w/8.;

      v2x = 3.-2.*r_2_stag-q2; 
      v2y = 3.-2.*r_2_stag+q2; 
      v2z = 1.+2.*r_2_stag+q2; 
      v2w = 1.+2.*r_2_stag-q2; 
      v2x = v2x/8.;
      v2y = v2y/8.;
      v2z = v2z/8.;
      v2w = v2w/8.;

      v3x = 3.-2.*r_3_c-q2; 
      v3y = 3.-2.*r_3_c+q2; 
      v3z = 1.+2.*r_3_c+q2; 
      v3w = 1.+2.*r_3_c-q2; 
      v3x = v3x/8.;
      v3y = v3y/8.;
      v3z = v3z/8.;
      v3w = v3w/8.;
      
      ws01 = v1x*v2x*v3x;
      ws11 = v1y*v2x*v3x;
      ws21 = v1z*v2x*v3x;
      ws31 = v1w*v2x*v3x;
      ws41 = v1x*v2y*v3x;
      ws51 = v1y*v2y*v3x;
      ws61 = v1z*v2y*v3x;
      ws71 = v1w*v2y*v3x;
      ws81 = v1x*v2z*v3x;
      ws91 = v1y*v2z*v3x;
      wsa1 = v1z*v2z*v3x;
      wsb1 = v1w*v2z*v3x;
      wsc1 = v1x*v2w*v3x;
      wsd1 = v1y*v2w*v3x;
      wse1 = v1z*v2w*v3x;
      wsf1 = v1w*v2w*v3x;
      
      ws02 = v1x*v2x*v3y;
      ws12 = v1y*v2x*v3y;
      ws22 = v1z*v2x*v3y;
      ws32 = v1w*v2x*v3y;
      ws42 = v1x*v2y*v3y;
      ws52 = v1y*v2y*v3y;
      ws62 = v1z*v2y*v3y;
      ws72 = v1w*v2y*v3y;
      ws82 = v1x*v2z*v3y;
      ws92 = v1y*v2z*v3y;
      wsa2 = v1z*v2z*v3y;
      wsb2 = v1w*v2z*v3y;
      wsc2 = v1x*v2w*v3y;
      wsd2 = v1y*v2w*v3y;
      wse2 = v1z*v2w*v3y;
      wsf2 = v1w*v2w*v3y;

      ws03 = v1x*v2x*v3z;
      ws13 = v1y*v2x*v3z;
      ws23 = v1z*v2x*v3z;
      ws33 = v1w*v2x*v3z;
      ws43 = v1x*v2y*v3z;
      ws53 = v1y*v2y*v3z;
      ws63 = v1z*v2y*v3z;
      ws73 = v1w*v2y*v3z;
      ws83 = v1x*v2z*v3z;
      ws93 = v1y*v2z*v3z;
      wsa3 = v1z*v2z*v3z;
      wsb3 = v1w*v2z*v3z;
      wsc3 = v1x*v2w*v3z;
      wsd3 = v1y*v2w*v3z;
      wse3 = v1z*v2w*v3z;
      wsf3 = v1w*v2w*v3z;

      ws04 = v1x*v2x*v3w;
      ws14 = v1y*v2x*v3w;
      ws24 = v1z*v2x*v3w;
      ws34 = v1w*v2x*v3w;
      ws44 = v1x*v2y*v3w;
      ws54 = v1y*v2y*v3w;
      ws64 = v1z*v2y*v3w;
      ws74 = v1w*v2y*v3w;
      ws84 = v1x*v2z*v3w;
      ws94 = v1y*v2z*v3w;
      wsa4 = v1z*v2z*v3w;
      wsb4 = v1w*v2z*v3w;
      wsc4 = v1x*v2w*v3w;
      wsd4 = v1y*v2w*v3w;
      wse4 = v1z*v2w*v3w;
      wsf4 = v1w*v2w*v3w;

      float my_F_2 = F_2[k2];
      f_2[ind_1x+N*ind_2x+N*N*ind_3x] += c*my_F_2*ws01;
      f_2[ind_1y+N*ind_2x+N*N*ind_3x] += c*my_F_2*ws11;
      f_2[ind_1z+N*ind_2x+N*N*ind_3x] += c*my_F_2*ws21;
      f_2[ind_1w+N*ind_2x+N*N*ind_3x] += c*my_F_2*ws31;
      f_2[ind_1x+N*ind_2y+N*N*ind_3x] += c*my_F_2*ws41;
      f_2[ind_1y+N*ind_2y+N*N*ind_3x] += c*my_F_2*ws51;
      f_2[ind_1z+N*ind_2y+N*N*ind_3x] += c*my_F_2*ws61;
      f_2[ind_1w+N*ind_2y+N*N*ind_3x] += c*my_F_2*ws71;
      f_2[ind_1x+N*ind_2z+N*N*ind_3x] += c*my_F_2*ws81;
      f_2[ind_1y+N*ind_2z+N*N*ind_3x] += c*my_F_2*ws91;
      f_2[ind_1z+N*ind_2z+N*N*ind_3x] += c*my_F_2*wsa1;
      f_2[ind_1w+N*ind_2z+N*N*ind_3x] += c*my_F_2*wsb1;
      f_2[ind_1x+N*ind_2w+N*N*ind_3x] += c*my_F_2*wsc1;
      f_2[ind_1y+N*ind_2w+N*N*ind_3x] += c*my_F_2*wsd1;
      f_2[ind_1z+N*ind_2w+N*N*ind_3x] += c*my_F_2*wse1;
      f_2[ind_1w+N*ind_2w+N*N*ind_3x] += c*my_F_2*wsf1;

      f_2[ind_1x+N*ind_2x+N*N*ind_3y] += c*my_F_2*ws02;
      f_2[ind_1y+N*ind_2x+N*N*ind_3y] += c*my_F_2*ws12;
      f_2[ind_1z+N*ind_2x+N*N*ind_3y] += c*my_F_2*ws22;
      f_2[ind_1w+N*ind_2x+N*N*ind_3y] += c*my_F_2*ws32;
      f_2[ind_1x+N*ind_2y+N*N*ind_3y] += c*my_F_2*ws42;
      f_2[ind_1y+N*ind_2y+N*N*ind_3y] += c*my_F_2*ws52;
      f_2[ind_1z+N*ind_2y+N*N*ind_3y] += c*my_F_2*ws62;
      f_2[ind_1w+N*ind_2y+N*N*ind_3y] += c*my_F_2*ws72;
      f_2[ind_1x+N*ind_2z+N*N*ind_3y] += c*my_F_2*ws82;
      f_2[ind_1y+N*ind_2z+N*N*ind_3y] += c*my_F_2*ws92;
      f_2[ind_1z+N*ind_2z+N*N*ind_3y] += c*my_F_2*wsa2;
      f_2[ind_1w+N*ind_2z+N*N*ind_3y] += c*my_F_2*wsb2;
      f_2[ind_1x+N*ind_2w+N*N*ind_3y] += c*my_F_2*wsc2;
      f_2[ind_1y+N*ind_2w+N*N*ind_3y] += c*my_F_2*wsd2;
      f_2[ind_1z+N*ind_2w+N*N*ind_3y] += c*my_F_2*wse2;
      f_2[ind_1w+N*ind_2w+N*N*ind_3y] += c*my_F_2*wsf2;

      f_2[ind_1x+N*ind_2x+N*N*ind_3z] += c*my_F_2*ws03;
      f_2[ind_1y+N*ind_2x+N*N*ind_3z] += c*my_F_2*ws13;
      f_2[ind_1z+N*ind_2x+N*N*ind_3z] += c*my_F_2*ws23;
      f_2[ind_1w+N*ind_2x+N*N*ind_3z] += c*my_F_2*ws33;
      f_2[ind_1x+N*ind_2y+N*N*ind_3z] += c*my_F_2*ws43;
      f_2[ind_1y+N*ind_2y+N*N*ind_3z] += c*my_F_2*ws53;
      f_2[ind_1z+N*ind_2y+N*N*ind_3z] += c*my_F_2*ws63;
      f_2[ind_1w+N*ind_2y+N*N*ind_3z] += c*my_F_2*ws73;
      f_2[ind_1x+N*ind_2z+N*N*ind_3z] += c*my_F_2*ws83;
      f_2[ind_1y+N*ind_2z+N*N*ind_3z] += c*my_F_2*ws93;
      f_2[ind_1z+N*ind_2z+N*N*ind_3z] += c*my_F_2*wsa3;
      f_2[ind_1w+N*ind_2z+N*N*ind_3z] += c*my_F_2*wsb3;
      f_2[ind_1x+N*ind_2w+N*N*ind_3z] += c*my_F_2*wsc3;
      f_2[ind_1y+N*ind_2w+N*N*ind_3z] += c*my_F_2*wsd3;
      f_2[ind_1z+N*ind_2w+N*N*ind_3z] += c*my_F_2*wse3;
      f_2[ind_1w+N*ind_2w+N*N*ind_3z] += c*my_F_2*wsf3;

      f_2[ind_1x+N*ind_2x+N*N*ind_3w] += c*my_F_2*ws04;
      f_2[ind_1y+N*ind_2x+N*N*ind_3w] += c*my_F_2*ws14;
      f_2[ind_1z+N*ind_2x+N*N*ind_3w] += c*my_F_2*ws24;
      f_2[ind_1w+N*ind_2x+N*N*ind_3w] += c*my_F_2*ws34;
      f_2[ind_1x+N*ind_2y+N*N*ind_3w] += c*my_F_2*ws44;
      f_2[ind_1y+N*ind_2y+N*N*ind_3w] += c*my_F_2*ws54;
      f_2[ind_1z+N*ind_2y+N*N*ind_3w] += c*my_F_2*ws64;
      f_2[ind_1w+N*ind_2y+N*N*ind_3w] += c*my_F_2*ws74;
      f_2[ind_1x+N*ind_2z+N*N*ind_3w] += c*my_F_2*ws84;
      f_2[ind_1y+N*ind_2z+N*N*ind_3w] += c*my_F_2*ws94;
      f_2[ind_1z+N*ind_2z+N*N*ind_3w] += c*my_F_2*wsa4;
      f_2[ind_1w+N*ind_2z+N*N*ind_3w] += c*my_F_2*wsb4;
      f_2[ind_1x+N*ind_2w+N*N*ind_3w] += c*my_F_2*wsc4;
      f_2[ind_1y+N*ind_2w+N*N*ind_3w] += c*my_F_2*wsd4;
      f_2[ind_1z+N*ind_2w+N*N*ind_3w] += c*my_F_2*wse4;
      f_2[ind_1w+N*ind_2w+N*N*ind_3w] += c*my_F_2*wsf4;

      k2 = nextn2[k2];
      }

      while(k3 != -1)
      {
      float s_1_stag = (X_1[k3]+0.5)/dx;
      float s_2_stag = (X_2[k3]+0.5)/dx;       
      float s_3_stag = (X_3[k3]+0.5)/dx;       
      int i_1_stag = convert_int(floor(s_1_stag));
      int i_2_stag = convert_int(floor(s_2_stag));
      int i_3_stag = convert_int(floor(s_3_stag));
      float r_1_stag = s_1_stag-i_1_stag;
      float r_2_stag = s_2_stag-i_2_stag;
      float r_3_stag = s_3_stag-i_3_stag;
      float s_1_c = (X_1[k3]+0.5-dx/2)/dx;
      float s_2_c = (X_2[k3]+0.5-dx/2)/dx;        
      float s_3_c = (X_3[k3]+0.5-dx/2)/dx;        
      int i_1_c = convert_int(floor(s_1_c));
      int i_2_c = convert_int(floor(s_2_c));
      int i_3_c = convert_int(floor(s_3_c));
      float r_1_c = s_1_c-i_1_c;
      float r_2_c = s_2_c-i_2_c;
      float r_3_c = s_3_c-i_3_c;
      
      ind_1x = (i_1_c-1 & N-1);
      ind_1y = (i_1_c & N-1);
      ind_1z = (i_1_c+1 & N-1);
      ind_1w = (i_1_c+2 & N-1);

      ind_2x = (i_2_c-1 & N-1);
      ind_2y = (i_2_c & N-1);
      ind_2z = (i_2_c+1 & N-1);
      ind_2w = (i_2_c+2 & N-1);

      ind_3x = (i_3_stag-1 & N-1);
      ind_3y = (i_3_stag & N-1);
      ind_3z = (i_3_stag+1 & N-1);
      ind_3w = (i_3_stag+2 & N-1);

      q1 = sqrt(1.+4.0*r_1_c*(1-r_1_c)); 
      q2 = sqrt(1.+4.0*r_2_c*(1-r_2_c));
      q3 = sqrt(1.+4.0*r_3_stag*(1-r_3_stag));

      v1x = 3.-2.*r_1_c-q1; 
      v1y = 3.-2.*r_1_c+q1; 
      v1z = 1.+2.*r_1_c+q1; 
      v1w = 1.+2.*r_1_c-q1; 
      v1x = v1x/8.;
      v1y = v1y/8.;
      v1z = v1z/8.;
      v1w = v1w/8.;

      v2x = 3.-2.*r_2_c-q2; 
      v2y = 3.-2.*r_2_c+q2; 
      v2z = 1.+2.*r_2_c+q2; 
      v2w = 1.+2.*r_2_c-q2; 
      v2x = v2x/8.;
      v2y = v2y/8.;
      v2z = v2z/8.;
      v2w = v2w/8.;

      v3x = 3.-2.*r_3_stag-q2; 
      v3y = 3.-2.*r_3_stag+q2; 
      v3z = 1.+2.*r_3_stag+q2; 
      v3w = 1.+2.*r_3_stag-q2; 
      v3x = v3x/8.;
      v3y = v3y/8.;
      v3z = v3z/8.;
      v3w = v3w/8.;
      
      ws01 = v1x*v2x*v3x;
      ws11 = v1y*v2x*v3x;
      ws21 = v1z*v2x*v3x;
      ws31 = v1w*v2x*v3x;
      ws41 = v1x*v2y*v3x;
      ws51 = v1y*v2y*v3x;
      ws61 = v1z*v2y*v3x;
      ws71 = v1w*v2y*v3x;
      ws81 = v1x*v2z*v3x;
      ws91 = v1y*v2z*v3x;
      wsa1 = v1z*v2z*v3x;
      wsb1 = v1w*v2z*v3x;
      wsc1 = v1x*v2w*v3x;
      wsd1 = v1y*v2w*v3x;
      wse1 = v1z*v2w*v3x;
      wsf1 = v1w*v2w*v3x;
      
      ws02 = v1x*v2x*v3y;
      ws12 = v1y*v2x*v3y;
      ws22 = v1z*v2x*v3y;
      ws32 = v1w*v2x*v3y;
      ws42 = v1x*v2y*v3y;
      ws52 = v1y*v2y*v3y;
      ws62 = v1z*v2y*v3y;
      ws72 = v1w*v2y*v3y;
      ws82 = v1x*v2z*v3y;
      ws92 = v1y*v2z*v3y;
      wsa2 = v1z*v2z*v3y;
      wsb2 = v1w*v2z*v3y;
      wsc2 = v1x*v2w*v3y;
      wsd2 = v1y*v2w*v3y;
      wse2 = v1z*v2w*v3y;
      wsf2 = v1w*v2w*v3y;

      ws03 = v1x*v2x*v3z;
      ws13 = v1y*v2x*v3z;
      ws23 = v1z*v2x*v3z;
      ws33 = v1w*v2x*v3z;
      ws43 = v1x*v2y*v3z;
      ws53 = v1y*v2y*v3z;
      ws63 = v1z*v2y*v3z;
      ws73 = v1w*v2y*v3z;
      ws83 = v1x*v2z*v3z;
      ws93 = v1y*v2z*v3z;
      wsa3 = v1z*v2z*v3z;
      wsb3 = v1w*v2z*v3z;
      wsc3 = v1x*v2w*v3z;
      wsd3 = v1y*v2w*v3z;
      wse3 = v1z*v2w*v3z;
      wsf3 = v1w*v2w*v3z;

      ws04 = v1x*v2x*v3w;
      ws14 = v1y*v2x*v3w;
      ws24 = v1z*v2x*v3w;
      ws34 = v1w*v2x*v3w;
      ws44 = v1x*v2y*v3w;
      ws54 = v1y*v2y*v3w;
      ws64 = v1z*v2y*v3w;
      ws74 = v1w*v2y*v3w;
      ws84 = v1x*v2z*v3w;
      ws94 = v1y*v2z*v3w;
      wsa4 = v1z*v2z*v3w;
      wsb4 = v1w*v2z*v3w;
      wsc4 = v1x*v2w*v3w;
      wsd4 = v1y*v2w*v3w;
      wse4 = v1z*v2w*v3w;
      wsf4 = v1w*v2w*v3w;

      float my_F_3 = F_3[k3];
      f_3[ind_1x+N*ind_2x+N*N*ind_3x] += c*my_F_3*ws01;
      f_3[ind_1y+N*ind_2x+N*N*ind_3x] += c*my_F_3*ws11;
      f_3[ind_1z+N*ind_2x+N*N*ind_3x] += c*my_F_3*ws21;
      f_3[ind_1w+N*ind_2x+N*N*ind_3x] += c*my_F_3*ws31;
      f_3[ind_1x+N*ind_2y+N*N*ind_3x] += c*my_F_3*ws41;
      f_3[ind_1y+N*ind_2y+N*N*ind_3x] += c*my_F_3*ws51;
      f_3[ind_1z+N*ind_2y+N*N*ind_3x] += c*my_F_3*ws61;
      f_3[ind_1w+N*ind_2y+N*N*ind_3x] += c*my_F_3*ws71;
      f_3[ind_1x+N*ind_2z+N*N*ind_3x] += c*my_F_3*ws81;
      f_3[ind_1y+N*ind_2z+N*N*ind_3x] += c*my_F_3*ws91;
      f_3[ind_1z+N*ind_2z+N*N*ind_3x] += c*my_F_3*wsa1;
      f_3[ind_1w+N*ind_2z+N*N*ind_3x] += c*my_F_3*wsb1;
      f_3[ind_1x+N*ind_2w+N*N*ind_3x] += c*my_F_3*wsc1;
      f_3[ind_1y+N*ind_2w+N*N*ind_3x] += c*my_F_3*wsd1;
      f_3[ind_1z+N*ind_2w+N*N*ind_3x] += c*my_F_3*wse1;
      f_3[ind_1w+N*ind_2w+N*N*ind_3x] += c*my_F_3*wsf1;

      f_3[ind_1x+N*ind_2x+N*N*ind_3y] += c*my_F_3*ws02;
      f_3[ind_1y+N*ind_2x+N*N*ind_3y] += c*my_F_3*ws12;
      f_3[ind_1z+N*ind_2x+N*N*ind_3y] += c*my_F_3*ws22;
      f_3[ind_1w+N*ind_2x+N*N*ind_3y] += c*my_F_3*ws32;
      f_3[ind_1x+N*ind_2y+N*N*ind_3y] += c*my_F_3*ws42;
      f_3[ind_1y+N*ind_2y+N*N*ind_3y] += c*my_F_3*ws52;
      f_3[ind_1z+N*ind_2y+N*N*ind_3y] += c*my_F_3*ws62;
      f_3[ind_1w+N*ind_2y+N*N*ind_3y] += c*my_F_3*ws72;
      f_3[ind_1x+N*ind_2z+N*N*ind_3y] += c*my_F_3*ws82;
      f_3[ind_1y+N*ind_2z+N*N*ind_3y] += c*my_F_3*ws92;
      f_3[ind_1z+N*ind_2z+N*N*ind_3y] += c*my_F_3*wsa2;
      f_3[ind_1w+N*ind_2z+N*N*ind_3y] += c*my_F_3*wsb2;
      f_3[ind_1x+N*ind_2w+N*N*ind_3y] += c*my_F_3*wsc2;
      f_3[ind_1y+N*ind_2w+N*N*ind_3y] += c*my_F_3*wsd2;
      f_3[ind_1z+N*ind_2w+N*N*ind_3y] += c*my_F_3*wse2;
      f_3[ind_1w+N*ind_2w+N*N*ind_3y] += c*my_F_3*wsf2;

      f_3[ind_1x+N*ind_2x+N*N*ind_3z] += c*my_F_3*ws03;
      f_3[ind_1y+N*ind_2x+N*N*ind_3z] += c*my_F_3*ws13;
      f_3[ind_1z+N*ind_2x+N*N*ind_3z] += c*my_F_3*ws23;
      f_3[ind_1w+N*ind_2x+N*N*ind_3z] += c*my_F_3*ws33;
      f_3[ind_1x+N*ind_2y+N*N*ind_3z] += c*my_F_3*ws43;
      f_3[ind_1y+N*ind_2y+N*N*ind_3z] += c*my_F_3*ws53;
      f_3[ind_1z+N*ind_2y+N*N*ind_3z] += c*my_F_3*ws63;
      f_3[ind_1w+N*ind_2y+N*N*ind_3z] += c*my_F_3*ws73;
      f_3[ind_1x+N*ind_2z+N*N*ind_3z] += c*my_F_3*ws83;
      f_3[ind_1y+N*ind_2z+N*N*ind_3z] += c*my_F_3*ws93;
      f_3[ind_1z+N*ind_2z+N*N*ind_3z] += c*my_F_3*wsa3;
      f_3[ind_1w+N*ind_2z+N*N*ind_3z] += c*my_F_3*wsb3;
      f_3[ind_1x+N*ind_2w+N*N*ind_3z] += c*my_F_3*wsc3;
      f_3[ind_1y+N*ind_2w+N*N*ind_3z] += c*my_F_3*wsd3;
      f_3[ind_1z+N*ind_2w+N*N*ind_3z] += c*my_F_3*wse3;
      f_3[ind_1w+N*ind_2w+N*N*ind_3z] += c*my_F_3*wsf3;

      f_3[ind_1x+N*ind_2x+N*N*ind_3w] += c*my_F_3*ws04;
      f_3[ind_1y+N*ind_2x+N*N*ind_3w] += c*my_F_3*ws14;
      f_3[ind_1z+N*ind_2x+N*N*ind_3w] += c*my_F_3*ws24;
      f_3[ind_1w+N*ind_2x+N*N*ind_3w] += c*my_F_3*ws34;
      f_3[ind_1x+N*ind_2y+N*N*ind_3w] += c*my_F_3*ws44;
      f_3[ind_1y+N*ind_2y+N*N*ind_3w] += c*my_F_3*ws54;
      f_3[ind_1z+N*ind_2y+N*N*ind_3w] += c*my_F_3*ws64;
      f_3[ind_1w+N*ind_2y+N*N*ind_3w] += c*my_F_3*ws74;
      f_3[ind_1x+N*ind_2z+N*N*ind_3w] += c*my_F_3*ws84;
      f_3[ind_1y+N*ind_2z+N*N*ind_3w] += c*my_F_3*ws94;
      f_3[ind_1z+N*ind_2z+N*N*ind_3w] += c*my_F_3*wsa4;
      f_3[ind_1w+N*ind_2z+N*N*ind_3w] += c*my_F_3*wsb4;
      f_3[ind_1x+N*ind_2w+N*N*ind_3w] += c*my_F_3*wsc4;
      f_3[ind_1y+N*ind_2w+N*N*ind_3w] += c*my_F_3*wsd4;
      f_3[ind_1z+N*ind_2w+N*N*ind_3w] += c*my_F_3*wse4;
      f_3[ind_1w+N*ind_2w+N*N*ind_3w] += c*my_F_3*wsf4;

      //k3 = nextn3[k3];
      }
    }
    __kernel void interp(__global const float *u_1, __global const float *u_2, __global const float *u_3, __global const float *X_1, __global const float *X_2, __global const float *X_3, const unsigned N, const float dx, __global float *U_1, __global float *U_2, __global float *U_3)
    {
      unsigned k = get_global_id(0);
      float s_1_stag = (X_1[k]+0.5)/dx;
      float s_2_stag = (X_2[k]+0.5)/dx;       
      float s_3_stag = (X_3[k]+0.5)/dx;       
      int i_1_stag = convert_int(floor(s_1_stag));
      int i_2_stag = convert_int(floor(s_2_stag));
      int i_3_stag = convert_int(floor(s_3_stag));
      float r_1_stag = s_1_stag-i_1_stag;
      float r_2_stag = s_2_stag-i_2_stag;
      float r_3_stag = s_3_stag-i_3_stag;
      float s_1_c = (X_1[k]+0.5-dx/2)/dx;
      float s_2_c = (X_2[k]+0.5-dx/2)/dx;        
      float s_3_c = (X_3[k]+0.5-dx/2)/dx;        
      int i_1_c = convert_int(floor(s_1_c));
      int i_2_c = convert_int(floor(s_2_c));
      int i_3_c = convert_int(floor(s_3_c));
      float r_1_c = s_1_c-i_1_c;
      float r_2_c = s_2_c-i_2_c;
      float r_3_c = s_3_c-i_3_c;

      int ind_1x, ind_1y, ind_1z, ind_1w;
      int ind_2x, ind_2y, ind_2z, ind_2w;
      int ind_3x, ind_3y, ind_3z, ind_3w;
      ind_1x = (i_1_stag-1 & N-1);
      ind_1y = (i_1_stag & N-1);
      ind_1z = (i_1_stag+1 & N-1);
      ind_1w = (i_1_stag+2 & N-1);

      ind_2x = (i_2_c-1 & N-1);
      ind_2y = (i_2_c & N-1);
      ind_2z = (i_2_c+1 & N-1);
      ind_2w = (i_2_c+2 & N-1);

      ind_3x = (i_3_c-1 & N-1);
      ind_3y = (i_3_c & N-1);
      ind_3z = (i_3_c+1 & N-1);
      ind_3w = (i_3_c+2 & N-1);

      float sb_tests01, sb_tests11, sb_tests21, sb_tests31;
      float sb_tests41, sb_tests51, sb_tests61, sb_tests71;
      float sb_tests81, sb_tests91, sb_testsa1, sb_testsb1;
      float sb_testsc1, sb_testsd1, sb_testse1, sb_testsf1;

      float sb_tests02, sb_tests12, sb_tests22, sb_tests32;
      float sb_tests42, sb_tests52, sb_tests62, sb_tests72;
      float sb_tests82, sb_tests92, sb_testsa2, sb_testsb2;
      float sb_testsc2, sb_testsd2, sb_testse2, sb_testsf2;

      float sb_tests03, sb_tests13, sb_tests23, sb_tests33;
      float sb_tests43, sb_tests53, sb_tests63, sb_tests73;
      float sb_tests83, sb_tests93, sb_testsa3, sb_testsb3;
      float sb_testsc3, sb_testsd3, sb_testse3, sb_testsf3;

      float sb_tests04, sb_tests14, sb_tests24, sb_tests34;
      float sb_tests44, sb_tests54, sb_tests64, sb_tests74;
      float sb_tests84, sb_tests94, sb_testsa4, sb_testsb4;
      float sb_testsc4, sb_testsd4, sb_testse4, sb_testsf4;

      sb_tests01 = u_1[ind_1x+N*ind_2x+N*N*ind_3x];
      sb_tests11 = u_1[ind_1y+N*ind_2x+N*N*ind_3x];
      sb_tests21 = u_1[ind_1z+N*ind_2x+N*N*ind_3x];
      sb_tests31 = u_1[ind_1w+N*ind_2x+N*N*ind_3x];
      sb_tests41 = u_1[ind_1x+N*ind_2y+N*N*ind_3x];
      sb_tests51 = u_1[ind_1y+N*ind_2y+N*N*ind_3x];
      sb_tests61 = u_1[ind_1z+N*ind_2y+N*N*ind_3x];
      sb_tests71 = u_1[ind_1w+N*ind_2y+N*N*ind_3x];
      sb_tests81 = u_1[ind_1x+N*ind_2z+N*N*ind_3x];
      sb_tests91 = u_1[ind_1y+N*ind_2z+N*N*ind_3x];
      sb_testsa1 = u_1[ind_1z+N*ind_2z+N*N*ind_3x];
      sb_testsb1 = u_1[ind_1w+N*ind_2z+N*N*ind_3x];
      sb_testsc1 = u_1[ind_1x+N*ind_2w+N*N*ind_3x];
      sb_testsd1 = u_1[ind_1y+N*ind_2w+N*N*ind_3x];
      sb_testse1 = u_1[ind_1z+N*ind_2w+N*N*ind_3x];
      sb_testsf1 = u_1[ind_1w+N*ind_2w+N*N*ind_3x];

      sb_tests02 = u_1[ind_1x+N*ind_2x+N*N*ind_3y];
      sb_tests12 = u_1[ind_1y+N*ind_2x+N*N*ind_3y];
      sb_tests22 = u_1[ind_1z+N*ind_2x+N*N*ind_3y];
      sb_tests32 = u_1[ind_1w+N*ind_2x+N*N*ind_3y];
      sb_tests42 = u_1[ind_1x+N*ind_2y+N*N*ind_3y];
      sb_tests52 = u_1[ind_1y+N*ind_2y+N*N*ind_3y];
      sb_tests62 = u_1[ind_1z+N*ind_2y+N*N*ind_3y];
      sb_tests72 = u_1[ind_1w+N*ind_2y+N*N*ind_3y];
      sb_tests82 = u_1[ind_1x+N*ind_2z+N*N*ind_3y];
      sb_tests92 = u_1[ind_1y+N*ind_2z+N*N*ind_3y];
      sb_testsa2 = u_1[ind_1z+N*ind_2z+N*N*ind_3y];
      sb_testsb2 = u_1[ind_1w+N*ind_2z+N*N*ind_3y];
      sb_testsc2 = u_1[ind_1x+N*ind_2w+N*N*ind_3y];
      sb_testsd2 = u_1[ind_1y+N*ind_2w+N*N*ind_3y];
      sb_testse2 = u_1[ind_1z+N*ind_2w+N*N*ind_3y];
      sb_testsf2 = u_1[ind_1w+N*ind_2w+N*N*ind_3y];

      sb_tests03 = u_1[ind_1x+N*ind_2x+N*N*ind_3z];
      sb_tests13 = u_1[ind_1y+N*ind_2x+N*N*ind_3z];
      sb_tests23 = u_1[ind_1z+N*ind_2x+N*N*ind_3z];
      sb_tests33 = u_1[ind_1w+N*ind_2x+N*N*ind_3z];
      sb_tests43 = u_1[ind_1x+N*ind_2y+N*N*ind_3z];
      sb_tests53 = u_1[ind_1y+N*ind_2y+N*N*ind_3z];
      sb_tests63 = u_1[ind_1z+N*ind_2y+N*N*ind_3z];
      sb_tests73 = u_1[ind_1w+N*ind_2y+N*N*ind_3z];
      sb_tests83 = u_1[ind_1x+N*ind_2z+N*N*ind_3z];
      sb_tests93 = u_1[ind_1y+N*ind_2z+N*N*ind_3z];
      sb_testsa3 = u_1[ind_1z+N*ind_2z+N*N*ind_3z];
      sb_testsb3 = u_1[ind_1w+N*ind_2z+N*N*ind_3z];
      sb_testsc3 = u_1[ind_1x+N*ind_2w+N*N*ind_3z];
      sb_testsd3 = u_1[ind_1y+N*ind_2w+N*N*ind_3z];
      sb_testse3 = u_1[ind_1z+N*ind_2w+N*N*ind_3z];
      sb_testsf3 = u_1[ind_1w+N*ind_2w+N*N*ind_3z];

      sb_tests04 = u_1[ind_1x+N*ind_2x+N*N*ind_3w];
      sb_tests14 = u_1[ind_1y+N*ind_2x+N*N*ind_3w];
      sb_tests24 = u_1[ind_1z+N*ind_2x+N*N*ind_3w];
      sb_tests34 = u_1[ind_1w+N*ind_2x+N*N*ind_3w];
      sb_tests44 = u_1[ind_1x+N*ind_2y+N*N*ind_3w];
      sb_tests54 = u_1[ind_1y+N*ind_2y+N*N*ind_3w];
      sb_tests64 = u_1[ind_1z+N*ind_2y+N*N*ind_3w];
      sb_tests74 = u_1[ind_1w+N*ind_2y+N*N*ind_3w];
      sb_tests84 = u_1[ind_1x+N*ind_2z+N*N*ind_3w];
      sb_tests94 = u_1[ind_1y+N*ind_2z+N*N*ind_3w];
      sb_testsa4 = u_1[ind_1z+N*ind_2z+N*N*ind_3w];
      sb_testsb4 = u_1[ind_1w+N*ind_2z+N*N*ind_3w];
      sb_testsc4 = u_1[ind_1x+N*ind_2w+N*N*ind_3w];
      sb_testsd4 = u_1[ind_1y+N*ind_2w+N*N*ind_3w];
      sb_testse4 = u_1[ind_1z+N*ind_2w+N*N*ind_3w];
      sb_testsf4 = u_1[ind_1w+N*ind_2w+N*N*ind_3w];

      float q1 = sqrt(1.+4.0*r_1_stag*(1-r_1_stag)); 
      float q2 = sqrt(1.+4.0*r_2_c*(1-r_2_c));
      float q3 = sqrt(1.+4.0*r_3_c*(1-r_3_c));

      float v1x, v1y, v1z, v1w;
      float v2x, v2y, v2z, v2w;
      float v3x, v3y, v3z, v3w;

      v1x = 3.-2.*r_1_stag-q1; 
      v1y = 3.-2.*r_1_stag+q1; 
      v1z = 1.+2.*r_1_stag+q1; 
      v1w = 1.+2.*r_1_stag-q1; 
      v1x = v1x/8.;
      v1y = v1y/8.;
      v1z = v1z/8.;
      v1w = v1w/8.;

      v2x = 3.-2.*r_2_c-q2; 
      v2y = 3.-2.*r_2_c+q2; 
      v2z = 1.+2.*r_2_c+q2; 
      v2w = 1.+2.*r_2_c-q2; 
      v2x = v2x/8.;
      v2y = v2y/8.;
      v2z = v2z/8.;
      v2w = v2w/8.;

      v3x = 3.-2.*r_3_c-q2; 
      v3y = 3.-2.*r_3_c+q2; 
      v3z = 1.+2.*r_3_c+q2; 
      v3w = 1.+2.*r_3_c-q2; 
      v3x = v3x/8.;
      v3y = v3y/8.;
      v3z = v3z/8.;
      v3w = v3w/8.;
      
      float ws01, ws11, ws21, ws31;
      float ws41, ws51, ws61, ws71;
      float ws81, ws91, wsa1, wsb1;
      float wsc1, wsd1, wse1, wsf1;

      float ws02, ws12, ws22, ws32;
      float ws42, ws52, ws62, ws72;
      float ws82, ws92, wsa2, wsb2;
      float wsc2, wsd2, wse2, wsf2;

      float ws03, ws13, ws23, ws33;
      float ws43, ws53, ws63, ws73;
      float ws83, ws93, wsa3, wsb3;
      float wsc3, wsd3, wse3, wsf3;

      float ws04, ws14, ws24, ws34;
      float ws44, ws54, ws64, ws74;
      float ws84, ws94, wsa4, wsb4;
      float wsc4, wsd4, wse4, wsf4;

      ws01 = v1x*v2x*v3x;
      ws11 = v1y*v2x*v3x;
      ws21 = v1z*v2x*v3x;
      ws31 = v1w*v2x*v3x;
      ws41 = v1x*v2y*v3x;
      ws51 = v1y*v2y*v3x;
      ws61 = v1z*v2y*v3x;
      ws71 = v1w*v2y*v3x;
      ws81 = v1x*v2z*v3x;
      ws91 = v1y*v2z*v3x;
      wsa1 = v1z*v2z*v3x;
      wsb1 = v1w*v2z*v3x;
      wsc1 = v1x*v2w*v3x;
      wsd1 = v1y*v2w*v3x;
      wse1 = v1z*v2w*v3x;
      wsf1 = v1w*v2w*v3x;
      
      ws02 = v1x*v2x*v3y;
      ws12 = v1y*v2x*v3y;
      ws22 = v1z*v2x*v3y;
      ws32 = v1w*v2x*v3y;
      ws42 = v1x*v2y*v3y;
      ws52 = v1y*v2y*v3y;
      ws62 = v1z*v2y*v3y;
      ws72 = v1w*v2y*v3y;
      ws82 = v1x*v2z*v3y;
      ws92 = v1y*v2z*v3y;
      wsa2 = v1z*v2z*v3y;
      wsb2 = v1w*v2z*v3y;
      wsc2 = v1x*v2w*v3y;
      wsd2 = v1y*v2w*v3y;
      wse2 = v1z*v2w*v3y;
      wsf2 = v1w*v2w*v3y;

      ws03 = v1x*v2x*v3z;
      ws13 = v1y*v2x*v3z;
      ws23 = v1z*v2x*v3z;
      ws33 = v1w*v2x*v3z;
      ws43 = v1x*v2y*v3z;
      ws53 = v1y*v2y*v3z;
      ws63 = v1z*v2y*v3z;
      ws73 = v1w*v2y*v3z;
      ws83 = v1x*v2z*v3z;
      ws93 = v1y*v2z*v3z;
      wsa3 = v1z*v2z*v3z;
      wsb3 = v1w*v2z*v3z;
      wsc3 = v1x*v2w*v3z;
      wsd3 = v1y*v2w*v3z;
      wse3 = v1z*v2w*v3z;
      wsf3 = v1w*v2w*v3z;

      ws04 = v1x*v2x*v3w;
      ws14 = v1y*v2x*v3w;
      ws24 = v1z*v2x*v3w;
      ws34 = v1w*v2x*v3w;
      ws44 = v1x*v2y*v3w;
      ws54 = v1y*v2y*v3w;
      ws64 = v1z*v2y*v3w;
      ws74 = v1w*v2y*v3w;
      ws84 = v1x*v2z*v3w;
      ws94 = v1y*v2z*v3w;
      wsa4 = v1z*v2z*v3w;
      wsb4 = v1w*v2z*v3w;
      wsc4 = v1x*v2w*v3w;
      wsd4 = v1y*v2w*v3w;
      wse4 = v1z*v2w*v3w;
      wsf4 = v1w*v2w*v3w;

      ws01 = ws01*sb_tests01;
      ws11 = ws11*sb_tests11;
      ws21 = ws21*sb_tests21;
      ws31 = ws31*sb_tests31;
      ws41 = ws41*sb_tests41;
      ws51 = ws51*sb_tests51;
      ws61 = ws61*sb_tests61;
      ws71 = ws71*sb_tests71;
      ws81 = ws81*sb_tests81;
      ws91 = ws91*sb_tests91;
      wsa1 = wsa1*sb_testsa1;
      wsb1 = wsb1*sb_testsb1;
      wsc1 = wsc1*sb_testsc1;
      wsd1 = wsd1*sb_testsd1;
      wse1 = wse1*sb_testse1;
      wsf1 = wsf1*sb_testsf1;

      ws02 = ws02*sb_tests02;
      ws12 = ws12*sb_tests12;
      ws22 = ws22*sb_tests22;
      ws32 = ws32*sb_tests32;
      ws42 = ws42*sb_tests42;
      ws52 = ws52*sb_tests52;
      ws62 = ws62*sb_tests62;
      ws72 = ws72*sb_tests72;
      ws82 = ws82*sb_tests82;
      ws92 = ws92*sb_tests92;
      wsa2 = wsa2*sb_testsa2;
      wsb2 = wsb2*sb_testsb2;
      wsc2 = wsc2*sb_testsc2;
      wsd2 = wsd2*sb_testsd2;
      wse2 = wse2*sb_testse2;
      wsf2 = wsf2*sb_testsf2;

      ws03 = ws03*sb_tests03;
      ws13 = ws13*sb_tests13;
      ws23 = ws23*sb_tests23;
      ws33 = ws33*sb_tests33;
      ws43 = ws43*sb_tests43;
      ws53 = ws53*sb_tests53;
      ws63 = ws63*sb_tests63;
      ws73 = ws73*sb_tests73;
      ws83 = ws83*sb_tests83;
      ws93 = ws93*sb_tests93;
      wsa3 = wsa3*sb_testsa3;
      wsb3 = wsb3*sb_testsb3;
      wsc3 = wsc3*sb_testsc3;
      wsd3 = wsd3*sb_testsd3;
      wse3 = wse3*sb_testse3;
      wsf3 = wsf3*sb_testsf3;

      ws04 = ws04*sb_tests04;
      ws14 = ws14*sb_tests14;
      ws24 = ws24*sb_tests24;
      ws34 = ws34*sb_tests34;
      ws44 = ws44*sb_tests44;
      ws54 = ws54*sb_tests54;
      ws64 = ws64*sb_tests64;
      ws74 = ws74*sb_tests74;
      ws84 = ws84*sb_tests84;
      ws94 = ws94*sb_tests94;
      wsa4 = wsa4*sb_testsa4;
      wsb4 = wsb4*sb_testsb4;
      wsc4 = wsc4*sb_testsc4;
      wsd4 = wsd4*sb_testsd4;
      wse4 = wse4*sb_testse4;
      wsf4 = wsf4*sb_testsf4;

      U_1[k] = ws01+ws11+ws21+ws31
	       +ws41+ws51+ws61+ws71
	       +ws81+ws91+wsa1+wsb1
	       +wsc1+wsd1+wse1+wsf1
               +ws02+ws12+ws22+ws32
	       +ws42+ws52+ws62+ws72
	       +ws82+ws92+wsa2+wsb2
	       +wsc2+wsd2+wse2+wsf2
               +ws03+ws13+ws23+ws33
	       +ws43+ws53+ws63+ws73
	       +ws83+ws93+wsa3+wsb3
	       +wsc3+wsd3+wse3+wsf3
               +ws04+ws14+ws24+ws34
	       +ws44+ws54+ws64+ws74
	       +ws84+ws94+wsa4+wsb4
	       +wsc4+wsd4+wse4+wsf4;

      ind_1x = (i_1_c-1 & N-1);
      ind_1y = (i_1_c & N-1);
      ind_1z = (i_1_c+1 & N-1);
      ind_1w = (i_1_c+2 & N-1);

      ind_2x = (i_2_stag-1 & N-1);
      ind_2y = (i_2_stag & N-1);
      ind_2z = (i_2_stag+1 & N-1);
      ind_2w = (i_2_stag+2 & N-1);

      ind_3x = (i_3_c-1 & N-1);
      ind_3y = (i_3_c & N-1);
      ind_3z = (i_3_c+1 & N-1);
      ind_3w = (i_3_c+2 & N-1);

      sb_tests01 = u_2[ind_1x+N*ind_2x+N*N*ind_3x];
      sb_tests11 = u_2[ind_1y+N*ind_2x+N*N*ind_3x];
      sb_tests21 = u_2[ind_1z+N*ind_2x+N*N*ind_3x];
      sb_tests31 = u_2[ind_1w+N*ind_2x+N*N*ind_3x];
      sb_tests41 = u_2[ind_1x+N*ind_2y+N*N*ind_3x];
      sb_tests51 = u_2[ind_1y+N*ind_2y+N*N*ind_3x];
      sb_tests61 = u_2[ind_1z+N*ind_2y+N*N*ind_3x];
      sb_tests71 = u_2[ind_1w+N*ind_2y+N*N*ind_3x];
      sb_tests81 = u_2[ind_1x+N*ind_2z+N*N*ind_3x];
      sb_tests91 = u_2[ind_1y+N*ind_2z+N*N*ind_3x];
      sb_testsa1 = u_2[ind_1z+N*ind_2z+N*N*ind_3x];
      sb_testsb1 = u_2[ind_1w+N*ind_2z+N*N*ind_3x];
      sb_testsc1 = u_2[ind_1x+N*ind_2w+N*N*ind_3x];
      sb_testsd1 = u_2[ind_1y+N*ind_2w+N*N*ind_3x];
      sb_testse1 = u_2[ind_1z+N*ind_2w+N*N*ind_3x];
      sb_testsf1 = u_2[ind_1w+N*ind_2w+N*N*ind_3x];

      sb_tests02 = u_2[ind_1x+N*ind_2x+N*N*ind_3y];
      sb_tests12 = u_2[ind_1y+N*ind_2x+N*N*ind_3y];
      sb_tests22 = u_2[ind_1z+N*ind_2x+N*N*ind_3y];
      sb_tests32 = u_2[ind_1w+N*ind_2x+N*N*ind_3y];
      sb_tests42 = u_2[ind_1x+N*ind_2y+N*N*ind_3y];
      sb_tests52 = u_2[ind_1y+N*ind_2y+N*N*ind_3y];
      sb_tests62 = u_2[ind_1z+N*ind_2y+N*N*ind_3y];
      sb_tests72 = u_2[ind_1w+N*ind_2y+N*N*ind_3y];
      sb_tests82 = u_2[ind_1x+N*ind_2z+N*N*ind_3y];
      sb_tests92 = u_2[ind_1y+N*ind_2z+N*N*ind_3y];
      sb_testsa2 = u_2[ind_1z+N*ind_2z+N*N*ind_3y];
      sb_testsb2 = u_2[ind_1w+N*ind_2z+N*N*ind_3y];
      sb_testsc2 = u_2[ind_1x+N*ind_2w+N*N*ind_3y];
      sb_testsd2 = u_2[ind_1y+N*ind_2w+N*N*ind_3y];
      sb_testse2 = u_2[ind_1z+N*ind_2w+N*N*ind_3y];
      sb_testsf2 = u_2[ind_1w+N*ind_2w+N*N*ind_3y];

      sb_tests03 = u_2[ind_1x+N*ind_2x+N*N*ind_3z];
      sb_tests13 = u_2[ind_1y+N*ind_2x+N*N*ind_3z];
      sb_tests23 = u_2[ind_1z+N*ind_2x+N*N*ind_3z];
      sb_tests33 = u_2[ind_1w+N*ind_2x+N*N*ind_3z];
      sb_tests43 = u_2[ind_1x+N*ind_2y+N*N*ind_3z];
      sb_tests53 = u_2[ind_1y+N*ind_2y+N*N*ind_3z];
      sb_tests63 = u_2[ind_1z+N*ind_2y+N*N*ind_3z];
      sb_tests73 = u_2[ind_1w+N*ind_2y+N*N*ind_3z];
      sb_tests83 = u_2[ind_1x+N*ind_2z+N*N*ind_3z];
      sb_tests93 = u_2[ind_1y+N*ind_2z+N*N*ind_3z];
      sb_testsa3 = u_2[ind_1z+N*ind_2z+N*N*ind_3z];
      sb_testsb3 = u_2[ind_1w+N*ind_2z+N*N*ind_3z];
      sb_testsc3 = u_2[ind_1x+N*ind_2w+N*N*ind_3z];
      sb_testsd3 = u_2[ind_1y+N*ind_2w+N*N*ind_3z];
      sb_testse3 = u_2[ind_1z+N*ind_2w+N*N*ind_3z];
      sb_testsf3 = u_2[ind_1w+N*ind_2w+N*N*ind_3z];

      sb_tests04 = u_2[ind_1x+N*ind_2x+N*N*ind_3w];
      sb_tests14 = u_2[ind_1y+N*ind_2x+N*N*ind_3w];
      sb_tests24 = u_2[ind_1z+N*ind_2x+N*N*ind_3w];
      sb_tests34 = u_2[ind_1w+N*ind_2x+N*N*ind_3w];
      sb_tests44 = u_2[ind_1x+N*ind_2y+N*N*ind_3w];
      sb_tests54 = u_2[ind_1y+N*ind_2y+N*N*ind_3w];
      sb_tests64 = u_2[ind_1z+N*ind_2y+N*N*ind_3w];
      sb_tests74 = u_2[ind_1w+N*ind_2y+N*N*ind_3w];
      sb_tests84 = u_2[ind_1x+N*ind_2z+N*N*ind_3w];
      sb_tests94 = u_2[ind_1y+N*ind_2z+N*N*ind_3w];
      sb_testsa4 = u_2[ind_1z+N*ind_2z+N*N*ind_3w];
      sb_testsb4 = u_2[ind_1w+N*ind_2z+N*N*ind_3w];
      sb_testsc4 = u_2[ind_1x+N*ind_2w+N*N*ind_3w];
      sb_testsd4 = u_2[ind_1y+N*ind_2w+N*N*ind_3w];
      sb_testse4 = u_2[ind_1z+N*ind_2w+N*N*ind_3w];
      sb_testsf4 = u_2[ind_1w+N*ind_2w+N*N*ind_3w];

      q1 = sqrt(1.+4.0*r_1_c*(1-r_1_c)); 
      q2 = sqrt(1.+4.0*r_2_stag*(1-r_2_stag));
      q3 = sqrt(1.+4.0*r_3_c*(1-r_3_c));

      v1x = 3.-2.*r_1_c-q1; 
      v1y = 3.-2.*r_1_c+q1; 
      v1z = 1.+2.*r_1_c+q1; 
      v1w = 1.+2.*r_1_c-q1; 
      v1x = v1x/8.;
      v1y = v1y/8.;
      v1z = v1z/8.;
      v1w = v1w/8.;

      v2x = 3.-2.*r_2_stag-q2; 
      v2y = 3.-2.*r_2_stag+q2; 
      v2z = 1.+2.*r_2_stag+q2; 
      v2w = 1.+2.*r_2_stag-q2; 
      v2x = v2x/8.;
      v2y = v2y/8.;
      v2z = v2z/8.;
      v2w = v2w/8.;

      v3x = 3.-2.*r_3_c-q2; 
      v3y = 3.-2.*r_3_c+q2; 
      v3z = 1.+2.*r_3_c+q2; 
      v3w = 1.+2.*r_3_c-q2; 
      v3x = v3x/8.;
      v3y = v3y/8.;
      v3z = v3z/8.;
      v3w = v3w/8.;
      
      ws01 = v1x*v2x*v3x;
      ws11 = v1y*v2x*v3x;
      ws21 = v1z*v2x*v3x;
      ws31 = v1w*v2x*v3x;
      ws41 = v1x*v2y*v3x;
      ws51 = v1y*v2y*v3x;
      ws61 = v1z*v2y*v3x;
      ws71 = v1w*v2y*v3x;
      ws81 = v1x*v2z*v3x;
      ws91 = v1y*v2z*v3x;
      wsa1 = v1z*v2z*v3x;
      wsb1 = v1w*v2z*v3x;
      wsc1 = v1x*v2w*v3x;
      wsd1 = v1y*v2w*v3x;
      wse1 = v1z*v2w*v3x;
      wsf1 = v1w*v2w*v3x;
      
      ws02 = v1x*v2x*v3y;
      ws12 = v1y*v2x*v3y;
      ws22 = v1z*v2x*v3y;
      ws32 = v1w*v2x*v3y;
      ws42 = v1x*v2y*v3y;
      ws52 = v1y*v2y*v3y;
      ws62 = v1z*v2y*v3y;
      ws72 = v1w*v2y*v3y;
      ws82 = v1x*v2z*v3y;
      ws92 = v1y*v2z*v3y;
      wsa2 = v1z*v2z*v3y;
      wsb2 = v1w*v2z*v3y;
      wsc2 = v1x*v2w*v3y;
      wsd2 = v1y*v2w*v3y;
      wse2 = v1z*v2w*v3y;
      wsf2 = v1w*v2w*v3y;

      ws03 = v1x*v2x*v3z;
      ws13 = v1y*v2x*v3z;
      ws23 = v1z*v2x*v3z;
      ws33 = v1w*v2x*v3z;
      ws43 = v1x*v2y*v3z;
      ws53 = v1y*v2y*v3z;
      ws63 = v1z*v2y*v3z;
      ws73 = v1w*v2y*v3z;
      ws83 = v1x*v2z*v3z;
      ws93 = v1y*v2z*v3z;
      wsa3 = v1z*v2z*v3z;
      wsb3 = v1w*v2z*v3z;
      wsc3 = v1x*v2w*v3z;
      wsd3 = v1y*v2w*v3z;
      wse3 = v1z*v2w*v3z;
      wsf3 = v1w*v2w*v3z;

      ws04 = v1x*v2x*v3w;
      ws14 = v1y*v2x*v3w;
      ws24 = v1z*v2x*v3w;
      ws34 = v1w*v2x*v3w;
      ws44 = v1x*v2y*v3w;
      ws54 = v1y*v2y*v3w;
      ws64 = v1z*v2y*v3w;
      ws74 = v1w*v2y*v3w;
      ws84 = v1x*v2z*v3w;
      ws94 = v1y*v2z*v3w;
      wsa4 = v1z*v2z*v3w;
      wsb4 = v1w*v2z*v3w;
      wsc4 = v1x*v2w*v3w;
      wsd4 = v1y*v2w*v3w;
      wse4 = v1z*v2w*v3w;
      wsf4 = v1w*v2w*v3w;

      ws01 = ws01*sb_tests01;
      ws11 = ws11*sb_tests11;
      ws21 = ws21*sb_tests21;
      ws31 = ws31*sb_tests31;
      ws41 = ws41*sb_tests41;
      ws51 = ws51*sb_tests51;
      ws61 = ws61*sb_tests61;
      ws71 = ws71*sb_tests71;
      ws81 = ws81*sb_tests81;
      ws91 = ws91*sb_tests91;
      wsa1 = wsa1*sb_testsa1;
      wsb1 = wsb1*sb_testsb1;
      wsc1 = wsc1*sb_testsc1;
      wsd1 = wsd1*sb_testsd1;
      wse1 = wse1*sb_testse1;
      wsf1 = wsf1*sb_testsf1;

      ws02 = ws02*sb_tests02;
      ws12 = ws12*sb_tests12;
      ws22 = ws22*sb_tests22;
      ws32 = ws32*sb_tests32;
      ws42 = ws42*sb_tests42;
      ws52 = ws52*sb_tests52;
      ws62 = ws62*sb_tests62;
      ws72 = ws72*sb_tests72;
      ws82 = ws82*sb_tests82;
      ws92 = ws92*sb_tests92;
      wsa2 = wsa2*sb_testsa2;
      wsb2 = wsb2*sb_testsb2;
      wsc2 = wsc2*sb_testsc2;
      wsd2 = wsd2*sb_testsd2;
      wse2 = wse2*sb_testse2;
      wsf2 = wsf2*sb_testsf2;

      ws03 = ws03*sb_tests03;
      ws13 = ws13*sb_tests13;
      ws23 = ws23*sb_tests23;
      ws33 = ws33*sb_tests33;
      ws43 = ws43*sb_tests43;
      ws53 = ws53*sb_tests53;
      ws63 = ws63*sb_tests63;
      ws73 = ws73*sb_tests73;
      ws83 = ws83*sb_tests83;
      ws93 = ws93*sb_tests93;
      wsa3 = wsa3*sb_testsa3;
      wsb3 = wsb3*sb_testsb3;
      wsc3 = wsc3*sb_testsc3;
      wsd3 = wsd3*sb_testsd3;
      wse3 = wse3*sb_testse3;
      wsf3 = wsf3*sb_testsf3;

      ws04 = ws04*sb_tests04;
      ws14 = ws14*sb_tests14;
      ws24 = ws24*sb_tests24;
      ws34 = ws34*sb_tests34;
      ws44 = ws44*sb_tests44;
      ws54 = ws54*sb_tests54;
      ws64 = ws64*sb_tests64;
      ws74 = ws74*sb_tests74;
      ws84 = ws84*sb_tests84;
      ws94 = ws94*sb_tests94;
      wsa4 = wsa4*sb_testsa4;
      wsb4 = wsb4*sb_testsb4;
      wsc4 = wsc4*sb_testsc4;
      wsd4 = wsd4*sb_testsd4;
      wse4 = wse4*sb_testse4;
      wsf4 = wsf4*sb_testsf4;

      U_2[k] = ws01+ws11+ws21+ws31
	       +ws41+ws51+ws61+ws71
	       +ws81+ws91+wsa1+wsb1
	       +wsc1+wsd1+wse1+wsf1
               +ws02+ws12+ws22+ws32
	       +ws42+ws52+ws62+ws72
	       +ws82+ws92+wsa2+wsb2
	       +wsc2+wsd2+wse2+wsf2
               +ws03+ws13+ws23+ws33
	       +ws43+ws53+ws63+ws73
	       +ws83+ws93+wsa3+wsb3
	       +wsc3+wsd3+wse3+wsf3
               +ws04+ws14+ws24+ws34
	       +ws44+ws54+ws64+ws74
	       +ws84+ws94+wsa4+wsb4
	       +wsc4+wsd4+wse4+wsf4;

      ind_1x = (i_1_c-1 & N-1);
      ind_1y = (i_1_c & N-1);
      ind_1z = (i_1_c+1 & N-1);
      ind_1w = (i_1_c+2 & N-1);

      ind_2x = (i_2_c-1 & N-1);
      ind_2y = (i_2_c & N-1);
      ind_2z = (i_2_c+1 & N-1);
      ind_2w = (i_2_c+2 & N-1);

      ind_3x = (i_3_stag-1 & N-1);
      ind_3y = (i_3_stag & N-1);
      ind_3z = (i_3_stag+1 & N-1);
      ind_3w = (i_3_stag+2 & N-1);

      sb_tests01 = u_3[ind_1x+N*ind_2x+N*N*ind_3x];
      sb_tests11 = u_3[ind_1y+N*ind_2x+N*N*ind_3x];
      sb_tests21 = u_3[ind_1z+N*ind_2x+N*N*ind_3x];
      sb_tests31 = u_3[ind_1w+N*ind_2x+N*N*ind_3x];
      sb_tests41 = u_3[ind_1x+N*ind_2y+N*N*ind_3x];
      sb_tests51 = u_3[ind_1y+N*ind_2y+N*N*ind_3x];
      sb_tests61 = u_3[ind_1z+N*ind_2y+N*N*ind_3x];
      sb_tests71 = u_3[ind_1w+N*ind_2y+N*N*ind_3x];
      sb_tests81 = u_3[ind_1x+N*ind_2z+N*N*ind_3x];
      sb_tests91 = u_3[ind_1y+N*ind_2z+N*N*ind_3x];
      sb_testsa1 = u_3[ind_1z+N*ind_2z+N*N*ind_3x];
      sb_testsb1 = u_3[ind_1w+N*ind_2z+N*N*ind_3x];
      sb_testsc1 = u_3[ind_1x+N*ind_2w+N*N*ind_3x];
      sb_testsd1 = u_3[ind_1y+N*ind_2w+N*N*ind_3x];
      sb_testse1 = u_3[ind_1z+N*ind_2w+N*N*ind_3x];
      sb_testsf1 = u_3[ind_1w+N*ind_2w+N*N*ind_3x];

      sb_tests02 = u_3[ind_1x+N*ind_2x+N*N*ind_3y];
      sb_tests12 = u_3[ind_1y+N*ind_2x+N*N*ind_3y];
      sb_tests22 = u_3[ind_1z+N*ind_2x+N*N*ind_3y];
      sb_tests32 = u_3[ind_1w+N*ind_2x+N*N*ind_3y];
      sb_tests42 = u_3[ind_1x+N*ind_2y+N*N*ind_3y];
      sb_tests52 = u_3[ind_1y+N*ind_2y+N*N*ind_3y];
      sb_tests62 = u_3[ind_1z+N*ind_2y+N*N*ind_3y];
      sb_tests72 = u_3[ind_1w+N*ind_2y+N*N*ind_3y];
      sb_tests82 = u_3[ind_1x+N*ind_2z+N*N*ind_3y];
      sb_tests92 = u_3[ind_1y+N*ind_2z+N*N*ind_3y];
      sb_testsa2 = u_3[ind_1z+N*ind_2z+N*N*ind_3y];
      sb_testsb2 = u_3[ind_1w+N*ind_2z+N*N*ind_3y];
      sb_testsc2 = u_3[ind_1x+N*ind_2w+N*N*ind_3y];
      sb_testsd2 = u_3[ind_1y+N*ind_2w+N*N*ind_3y];
      sb_testse2 = u_3[ind_1z+N*ind_2w+N*N*ind_3y];
      sb_testsf2 = u_3[ind_1w+N*ind_2w+N*N*ind_3y];

      sb_tests03 = u_3[ind_1x+N*ind_2x+N*N*ind_3z];
      sb_tests13 = u_3[ind_1y+N*ind_2x+N*N*ind_3z];
      sb_tests23 = u_3[ind_1z+N*ind_2x+N*N*ind_3z];
      sb_tests33 = u_3[ind_1w+N*ind_2x+N*N*ind_3z];
      sb_tests43 = u_3[ind_1x+N*ind_2y+N*N*ind_3z];
      sb_tests53 = u_3[ind_1y+N*ind_2y+N*N*ind_3z];
      sb_tests63 = u_3[ind_1z+N*ind_2y+N*N*ind_3z];
      sb_tests73 = u_3[ind_1w+N*ind_2y+N*N*ind_3z];
      sb_tests83 = u_3[ind_1x+N*ind_2z+N*N*ind_3z];
      sb_tests93 = u_3[ind_1y+N*ind_2z+N*N*ind_3z];
      sb_testsa3 = u_3[ind_1z+N*ind_2z+N*N*ind_3z];
      sb_testsb3 = u_3[ind_1w+N*ind_2z+N*N*ind_3z];
      sb_testsc3 = u_3[ind_1x+N*ind_2w+N*N*ind_3z];
      sb_testsd3 = u_3[ind_1y+N*ind_2w+N*N*ind_3z];
      sb_testse3 = u_3[ind_1z+N*ind_2w+N*N*ind_3z];
      sb_testsf3 = u_3[ind_1w+N*ind_2w+N*N*ind_3z];

      sb_tests04 = u_3[ind_1x+N*ind_2x+N*N*ind_3w];
      sb_tests14 = u_3[ind_1y+N*ind_2x+N*N*ind_3w];
      sb_tests24 = u_3[ind_1z+N*ind_2x+N*N*ind_3w];
      sb_tests34 = u_3[ind_1w+N*ind_2x+N*N*ind_3w];
      sb_tests44 = u_3[ind_1x+N*ind_2y+N*N*ind_3w];
      sb_tests54 = u_3[ind_1y+N*ind_2y+N*N*ind_3w];
      sb_tests64 = u_3[ind_1z+N*ind_2y+N*N*ind_3w];
      sb_tests74 = u_3[ind_1w+N*ind_2y+N*N*ind_3w];
      sb_tests84 = u_3[ind_1x+N*ind_2z+N*N*ind_3w];
      sb_tests94 = u_3[ind_1y+N*ind_2z+N*N*ind_3w];
      sb_testsa4 = u_3[ind_1z+N*ind_2z+N*N*ind_3w];
      sb_testsb4 = u_3[ind_1w+N*ind_2z+N*N*ind_3w];
      sb_testsc4 = u_3[ind_1x+N*ind_2w+N*N*ind_3w];
      sb_testsd4 = u_3[ind_1y+N*ind_2w+N*N*ind_3w];
      sb_testse4 = u_3[ind_1z+N*ind_2w+N*N*ind_3w];
      sb_testsf4 = u_3[ind_1w+N*ind_2w+N*N*ind_3w];

      q1 = sqrt(1.+4.0*r_1_c*(1-r_1_c)); 
      q2 = sqrt(1.+4.0*r_2_c*(1-r_2_c));
      q3 = sqrt(1.+4.0*r_3_stag*(1-r_3_stag));

      v1x = 3.-2.*r_1_c-q1; 
      v1y = 3.-2.*r_1_c+q1; 
      v1z = 1.+2.*r_1_c+q1; 
      v1w = 1.+2.*r_1_c-q1; 
      v1x = v1x/8.;
      v1y = v1y/8.;
      v1z = v1z/8.;
      v1w = v1w/8.;

      v2x = 3.-2.*r_2_c-q2; 
      v2y = 3.-2.*r_2_c+q2; 
      v2z = 1.+2.*r_2_c+q2; 
      v2w = 1.+2.*r_2_c-q2; 
      v2x = v2x/8.;
      v2y = v2y/8.;
      v2z = v2z/8.;
      v2w = v2w/8.;

      v3x = 3.-2.*r_3_stag-q2; 
      v3y = 3.-2.*r_3_stag+q2; 
      v3z = 1.+2.*r_3_stag+q2; 
      v3w = 1.+2.*r_3_stag-q2; 
      v3x = v3x/8.;
      v3y = v3y/8.;
      v3z = v3z/8.;
      v3w = v3w/8.;
      
      ws01 = v1x*v2x*v3x;
      ws11 = v1y*v2x*v3x;
      ws21 = v1z*v2x*v3x;
      ws31 = v1w*v2x*v3x;
      ws41 = v1x*v2y*v3x;
      ws51 = v1y*v2y*v3x;
      ws61 = v1z*v2y*v3x;
      ws71 = v1w*v2y*v3x;
      ws81 = v1x*v2z*v3x;
      ws91 = v1y*v2z*v3x;
      wsa1 = v1z*v2z*v3x;
      wsb1 = v1w*v2z*v3x;
      wsc1 = v1x*v2w*v3x;
      wsd1 = v1y*v2w*v3x;
      wse1 = v1z*v2w*v3x;
      wsf1 = v1w*v2w*v3x;
      
      ws02 = v1x*v2x*v3y;
      ws12 = v1y*v2x*v3y;
      ws22 = v1z*v2x*v3y;
      ws32 = v1w*v2x*v3y;
      ws42 = v1x*v2y*v3y;
      ws52 = v1y*v2y*v3y;
      ws62 = v1z*v2y*v3y;
      ws72 = v1w*v2y*v3y;
      ws82 = v1x*v2z*v3y;
      ws92 = v1y*v2z*v3y;
      wsa2 = v1z*v2z*v3y;
      wsb2 = v1w*v2z*v3y;
      wsc2 = v1x*v2w*v3y;
      wsd2 = v1y*v2w*v3y;
      wse2 = v1z*v2w*v3y;
      wsf2 = v1w*v2w*v3y;

      ws03 = v1x*v2x*v3z;
      ws13 = v1y*v2x*v3z;
      ws23 = v1z*v2x*v3z;
      ws33 = v1w*v2x*v3z;
      ws43 = v1x*v2y*v3z;
      ws53 = v1y*v2y*v3z;
      ws63 = v1z*v2y*v3z;
      ws73 = v1w*v2y*v3z;
      ws83 = v1x*v2z*v3z;
      ws93 = v1y*v2z*v3z;
      wsa3 = v1z*v2z*v3z;
      wsb3 = v1w*v2z*v3z;
      wsc3 = v1x*v2w*v3z;
      wsd3 = v1y*v2w*v3z;
      wse3 = v1z*v2w*v3z;
      wsf3 = v1w*v2w*v3z;

      ws04 = v1x*v2x*v3w;
      ws14 = v1y*v2x*v3w;
      ws24 = v1z*v2x*v3w;
      ws34 = v1w*v2x*v3w;
      ws44 = v1x*v2y*v3w;
      ws54 = v1y*v2y*v3w;
      ws64 = v1z*v2y*v3w;
      ws74 = v1w*v2y*v3w;
      ws84 = v1x*v2z*v3w;
      ws94 = v1y*v2z*v3w;
      wsa4 = v1z*v2z*v3w;
      wsb4 = v1w*v2z*v3w;
      wsc4 = v1x*v2w*v3w;
      wsd4 = v1y*v2w*v3w;
      wse4 = v1z*v2w*v3w;
      wsf4 = v1w*v2w*v3w;

      ws01 = ws01*sb_tests01;
      ws11 = ws11*sb_tests11;
      ws21 = ws21*sb_tests21;
      ws31 = ws31*sb_tests31;
      ws41 = ws41*sb_tests41;
      ws51 = ws51*sb_tests51;
      ws61 = ws61*sb_tests61;
      ws71 = ws71*sb_tests71;
      ws81 = ws81*sb_tests81;
      ws91 = ws91*sb_tests91;
      wsa1 = wsa1*sb_testsa1;
      wsb1 = wsb1*sb_testsb1;
      wsc1 = wsc1*sb_testsc1;
      wsd1 = wsd1*sb_testsd1;
      wse1 = wse1*sb_testse1;
      wsf1 = wsf1*sb_testsf1;

      ws02 = ws02*sb_tests02;
      ws12 = ws12*sb_tests12;
      ws22 = ws22*sb_tests22;
      ws32 = ws32*sb_tests32;
      ws42 = ws42*sb_tests42;
      ws52 = ws52*sb_tests52;
      ws62 = ws62*sb_tests62;
      ws72 = ws72*sb_tests72;
      ws82 = ws82*sb_tests82;
      ws92 = ws92*sb_tests92;
      wsa2 = wsa2*sb_testsa2;
      wsb2 = wsb2*sb_testsb2;
      wsc2 = wsc2*sb_testsc2;
      wsd2 = wsd2*sb_testsd2;
      wse2 = wse2*sb_testse2;
      wsf2 = wsf2*sb_testsf2;

      ws03 = ws03*sb_tests03;
      ws13 = ws13*sb_tests13;
      ws23 = ws23*sb_tests23;
      ws33 = ws33*sb_tests33;
      ws43 = ws43*sb_tests43;
      ws53 = ws53*sb_tests53;
      ws63 = ws63*sb_tests63;
      ws73 = ws73*sb_tests73;
      ws83 = ws83*sb_tests83;
      ws93 = ws93*sb_tests93;
      wsa3 = wsa3*sb_testsa3;
      wsb3 = wsb3*sb_testsb3;
      wsc3 = wsc3*sb_testsc3;
      wsd3 = wsd3*sb_testsd3;
      wse3 = wse3*sb_testse3;
      wsf3 = wsf3*sb_testsf3;

      ws04 = ws04*sb_tests04;
      ws14 = ws14*sb_tests14;
      ws24 = ws24*sb_tests24;
      ws34 = ws34*sb_tests34;
      ws44 = ws44*sb_tests44;
      ws54 = ws54*sb_tests54;
      ws64 = ws64*sb_tests64;
      ws74 = ws74*sb_tests74;
      ws84 = ws84*sb_tests84;
      ws94 = ws94*sb_tests94;
      wsa4 = wsa4*sb_testsa4;
      wsb4 = wsb4*sb_testsb4;
      wsc4 = wsc4*sb_testsc4;
      wsd4 = wsd4*sb_testsd4;
      wse4 = wse4*sb_testse4;
      wsf4 = wsf4*sb_testsf4;

      U_3[k] = ws01+ws11+ws21+ws31
	       +ws41+ws51+ws61+ws71
	       +ws81+ws91+wsa1+wsb1
	       +wsc1+wsd1+wse1+wsf1
               +ws02+ws12+ws22+ws32
	       +ws42+ws52+ws62+ws72
	       +ws82+ws92+wsa2+wsb2
	       +wsc2+wsd2+wse2+wsf2
               +ws03+ws13+ws23+ws33
	       +ws43+ws53+ws63+ws73
	       +ws83+ws93+wsa3+wsb3
	       +wsc3+wsd3+wse3+wsf3
               +ws04+ws14+ws24+ws34
	       +ws44+ws54+ws64+ws74
	       +ws84+ws94+wsa4+wsb4
	       +wsc4+wsd4+wse4+wsf4;
    }
    __kernel void spread_save(__global const int *firstn1, __global const int *nextn1, __global const int *firstn2, __global const int *nextn2, __global const float *F_1, __global const float *F_2, __global const float *X_1, __global const float *X_2, const unsigned N, const float dx, const float dtheta, const unsigned shift0, const unsigned shift1, __global float *f_1, __global float *f_2)
    {
      unsigned i = 4*get_global_id(0)+shift0;
      unsigned j = 4*get_global_id(1)+shift1;

      int k1 = firstn1[i+N*j]; 
      int k2 = firstn2[i+N*j]; 
      if (k1 == -1 && k2 == -1){
      return;
      }

      float c = dtheta/(dx*dx*dx);
      int ind_1x, ind_1y, ind_1z, ind_1w;
      int ind_2x, ind_2y, ind_2z, ind_2w;

      float f_10, f_11, f_12, f_13;
      float f_14, f_15, f_16, f_17;
      float f_18, f_19, f_1a, f_1b;
      float f_1c, f_1d, f_1e, f_1f;

      float f_20, f_21, f_22, f_23;
      float f_24, f_25, f_26, f_27;
      float f_28, f_29, f_2a, f_2b;
      float f_2c, f_2d, f_2e, f_2f;

      float v1x, v1y, v1z, v1w;
      float v2x, v2y, v2z, v2w;

      float q1, q2;

      float ws0, ws1, ws2, ws3;
      float ws4, ws5, ws6, ws7;
      float ws8, ws9, wsa, wsb;
      float wsc, wsd, wse, wsf;

      if (k1 !=-1){
      ind_1x = (i-1 & N-1);
      ind_1y = (i & N-1);
      ind_1z = (i+1 & N-1);
      ind_1w = (i+2 & N-1);

      ind_2x = (j-1 & N-1);
      ind_2y = (j & N-1);
      ind_2z = (j+1 & N-1);
      ind_2w = (j+2 & N-1);

      f_10=f_1[ind_1x+N*ind_2x];
      f_11=f_1[ind_1y+N*ind_2x];
      f_12=f_1[ind_1z+N*ind_2x];
      f_13=f_1[ind_1w+N*ind_2x];
      f_14=f_1[ind_1x+N*ind_2y];
      f_15=f_1[ind_1y+N*ind_2y];
      f_16=f_1[ind_1z+N*ind_2y];
      f_17=f_1[ind_1w+N*ind_2y];
      f_18=f_1[ind_1x+N*ind_2z];
      f_19=f_1[ind_1y+N*ind_2z];
      f_1a=f_1[ind_1z+N*ind_2z];
      f_1b=f_1[ind_1w+N*ind_2z];
      f_1c=f_1[ind_1x+N*ind_2w];
      f_1d=f_1[ind_1y+N*ind_2w];
      f_1e=f_1[ind_1z+N*ind_2w];
      f_1f=f_1[ind_1w+N*ind_2w];
      }

      while(k1 != -1)
      {
      float s_1_stag = (X_1[k1]+0.5)/dx;
      int i_1_stag = convert_int(floor(s_1_stag));
      float r_1_stag = s_1_stag-i_1_stag;
      float s_2_c = (X_2[k1]+0.5-dx/2)/dx;        
      int i_2_c = convert_int(floor(s_2_c));
      float r_2_c = s_2_c-i_2_c;
      
      q1 = sqrt(1.+4.0*r_1_stag*(1-r_1_stag)); 
      q2 = sqrt(1.+4.0*r_2_c*(1-r_2_c));

      v1x = 3.-2.*r_1_stag-q1; 
      v1y = 3.-2.*r_1_stag+q1; 
      v1z = 1.+2.*r_1_stag+q1; 
      v1w = 1.+2.*r_1_stag-q1; 
      v1x = v1x/8.;
      v1y = v1y/8.;
      v1z = v1z/8.;
      v1w = v1w/8.;

      v2x = 3.-2.*r_2_c-q2; 
      v2y = 3.-2.*r_2_c+q2; 
      v2z = 1.+2.*r_2_c+q2; 
      v2w = 1.+2.*r_2_c-q2; 
      v2x = v2x/8.;
      v2y = v2y/8.;
      v2z = v2z/8.;
      v2w = v2w/8.;
      
      ws0 = v1x*v2x;
      ws1 = v1y*v2x;
      ws2 = v1z*v2x;
      ws3 = v1w*v2x;
      ws4 = v1x*v2y;
      ws5 = v1y*v2y;
      ws6 = v1z*v2y;
      ws7 = v1w*v2y;
      ws8 = v1x*v2z;
      ws9 = v1y*v2z;
      wsa = v1z*v2z;
      wsb = v1w*v2z;
      wsc = v1x*v2w;
      wsd = v1y*v2w;
      wse = v1z*v2w;
      wsf = v1w*v2w;
      
      float my_F_1 = F_1[k1];
      f_10 += c*my_F_1*ws0;
      f_11 += c*my_F_1*ws1;
      f_12 += c*my_F_1*ws2;
      f_13 += c*my_F_1*ws3;
      f_14 += c*my_F_1*ws4;
      f_15 += c*my_F_1*ws5;
      f_16 += c*my_F_1*ws6;
      f_17 += c*my_F_1*ws7;
      f_18 += c*my_F_1*ws8;
      f_19 += c*my_F_1*ws9;
      f_1a += c*my_F_1*wsa;
      f_1b += c*my_F_1*wsb;
      f_1c += c*my_F_1*wsc;
      f_1d += c*my_F_1*wsd;
      f_1e += c*my_F_1*wse;
      f_1f += c*my_F_1*wsf;

      k1 = nextn1[k1];
      if (k1 == -1){
      f_1[ind_1x+N*ind_2x]=f_10;
      f_1[ind_1y+N*ind_2x]=f_11;
      f_1[ind_1z+N*ind_2x]=f_12;
      f_1[ind_1w+N*ind_2x]=f_13;
      f_1[ind_1x+N*ind_2y]=f_14;
      f_1[ind_1y+N*ind_2y]=f_15;
      f_1[ind_1z+N*ind_2y]=f_16;
      f_1[ind_1w+N*ind_2y]=f_17;
      f_1[ind_1x+N*ind_2z]=f_18;
      f_1[ind_1y+N*ind_2z]=f_19;
      f_1[ind_1z+N*ind_2z]=f_1a;
      f_1[ind_1w+N*ind_2z]=f_1b;
      f_1[ind_1x+N*ind_2w]=f_1c;
      f_1[ind_1y+N*ind_2w]=f_1d;
      f_1[ind_1z+N*ind_2w]=f_1e;
      f_1[ind_1w+N*ind_2w]=f_1f;
      }
      }

      if (k2 !=-1){
      ind_1x = (i-1 & N-1);
      ind_1y = (i & N-1);
      ind_1z = (i+1 & N-1);
      ind_1w = (i+2 & N-1);

      ind_2x = (j-1 & N-1);
      ind_2y = (j & N-1);
      ind_2z = (j+1 & N-1);
      ind_2w = (j+2 & N-1);

      f_20=f_2[ind_1x+N*ind_2x];
      f_21=f_2[ind_1y+N*ind_2x];
      f_22=f_2[ind_1z+N*ind_2x];
      f_23=f_2[ind_1w+N*ind_2x];
      f_24=f_2[ind_1x+N*ind_2y];
      f_25=f_2[ind_1y+N*ind_2y];
      f_26=f_2[ind_1z+N*ind_2y];
      f_27=f_2[ind_1w+N*ind_2y];
      f_28=f_2[ind_1x+N*ind_2z];
      f_29=f_2[ind_1y+N*ind_2z];
      f_2a=f_2[ind_1z+N*ind_2z];
      f_2b=f_2[ind_1w+N*ind_2z];
      f_2c=f_2[ind_1x+N*ind_2w];
      f_2d=f_2[ind_1y+N*ind_2w];
      f_2e=f_2[ind_1z+N*ind_2w];
      f_2f=f_2[ind_1w+N*ind_2w];
      }
      
      while (k2 != -1){
      float s_1_c = (X_1[k2]+0.5-dx/2)/dx;
      int i_1_c = convert_int(floor(s_1_c));
      float r_1_c = s_1_c-i_1_c;
      float s_2_stag = (X_2[k2]+0.5)/dx;        
      int i_2_stag = convert_int(floor(s_2_stag));
      float r_2_stag = s_2_stag-i_2_stag;

      q1 = sqrt(1.+4.0*r_1_c*(1-r_1_c)); 
      q2 = sqrt(1.+4.0*r_2_stag*(1-r_2_stag));

      v1x = 3.-2.*r_1_c-q1; 
      v1y = 3.-2.*r_1_c+q1; 
      v1z = 1.+2.*r_1_c+q1; 
      v1w = 1.+2.*r_1_c-q1; 
      v1x = v1x/8.;
      v1y = v1y/8.;
      v1z = v1z/8.;
      v1w = v1w/8.;

      v2x = 3.-2.*r_2_stag-q2; 
      v2y = 3.-2.*r_2_stag+q2; 
      v2z = 1.+2.*r_2_stag+q2; 
      v2w = 1.+2.*r_2_stag-q2; 
      v2x = v2x/8.;
      v2y = v2y/8.;
      v2z = v2z/8.;
      v2w = v2w/8.;
      
      ws0 = v1x*v2x;
      ws1 = v1y*v2x;
      ws2 = v1z*v2x;
      ws3 = v1w*v2x;
      ws4 = v1x*v2y;
      ws5 = v1y*v2y;
      ws6 = v1z*v2y;
      ws7 = v1w*v2y;
      ws8 = v1x*v2z;
      ws9 = v1y*v2z;
      wsa = v1z*v2z;
      wsb = v1w*v2z;
      wsc = v1x*v2w;
      wsd = v1y*v2w;
      wse = v1z*v2w;
      wsf = v1w*v2w;
      
      float my_F_2 = F_2[k2];
      f_20 += c*my_F_2*ws0;
      f_21 += c*my_F_2*ws1;
      f_22 += c*my_F_2*ws2;
      f_23 += c*my_F_2*ws3;
      f_24 += c*my_F_2*ws4;
      f_25 += c*my_F_2*ws5;
      f_26 += c*my_F_2*ws6;
      f_27 += c*my_F_2*ws7;
      f_28 += c*my_F_2*ws8;
      f_29 += c*my_F_2*ws9;
      f_2a += c*my_F_2*wsa;
      f_2b += c*my_F_2*wsb;
      f_2c += c*my_F_2*wsc;
      f_2d += c*my_F_2*wsd;
      f_2e += c*my_F_2*wse;
      f_2f += c*my_F_2*wsf;

      k2 = nextn2[k2];
      if (k2 == -1){
      f_2[ind_1x+N*ind_2x]=f_20;
      f_2[ind_1y+N*ind_2x]=f_21;
      f_2[ind_1z+N*ind_2x]=f_22;
      f_2[ind_1w+N*ind_2x]=f_23;
      f_2[ind_1x+N*ind_2y]=f_24;
      f_2[ind_1y+N*ind_2y]=f_25;
      f_2[ind_1z+N*ind_2y]=f_26;
      f_2[ind_1w+N*ind_2y]=f_27;
      f_2[ind_1x+N*ind_2z]=f_28;
      f_2[ind_1y+N*ind_2z]=f_29;
      f_2[ind_1z+N*ind_2z]=f_2a;
      f_2[ind_1w+N*ind_2z]=f_2b;
      f_2[ind_1x+N*ind_2w]=f_2c;
      f_2[ind_1y+N*ind_2w]=f_2d;
      f_2[ind_1z+N*ind_2w]=f_2e;
      f_2[ind_1w+N*ind_2w]=f_2f;
      }
      }
    }
    __kernel void updaten(__global int *firstn1, __global int *nextn1, __global int *firstn2, __global int *nextn2, __global const float *X_1, __global const float *X_2, const unsigned N, const float dx, const unsigned shift0, const unsigned shift1)
    {
      unsigned i = 4*get_global_id(0)+shift0;
      unsigned j = 4*get_global_id(1)+shift1;

      int k1 = firstn1[i+N*j]; 
      int k2 = firstn2[i+N*j]; 
      if (k1 == -1 && k2 == -1){
      return;
      }

      uint ind_1, ind_2;
      int nextnold;

      int nprev1 = -1;
      while(k1 != -1)
      {
      float s_1_stag = (X_1[k1]+0.5)/dx;
      int i_1_stag = convert_int(floor(s_1_stag));
      float s_2_c = (X_2[k1]+0.5-dx/2)/dx;        
      int i_2_c = convert_int(floor(s_2_c));
      
      ind_1 = (i_1_stag & N-1);
      ind_2 = (i_2_c & N-1);

      if (ind_1 == i && ind_2 == j){
      nprev1 = k1;
      k1 = nextn1[k1];
      }
      else {
      nextnold = nextn1[k1];
      nextn1[k1] = firstn1[ind_1+N*ind_2];
      firstn1[ind_1+N*ind_2] = k1;
      k1 = nextnold;
          if (nprev1 == -1) {
          firstn1[i+N*j] = nextnold;
          }
          else {
          nextn1[nprev1] = nextnold;
          }
      }
      }

      int nprev2 = -1;
      while(k2 != -1)
      {
      float s_2_stag = (X_2[k2]+0.5)/dx;        
      int i_2_stag = convert_int(floor(s_2_stag));
      float s_1_c = (X_1[k2]+0.5-dx/2)/dx;
      int i_1_c = convert_int(floor(s_1_c));
      
      ind_1 = (i_1_c & N-1);
      ind_2 = (i_2_stag & N-1);

      if (ind_1 == i && ind_2 == j){
      nprev2 = k2;
      k2 = nextn2[k2];
      }
      else {
      nextnold = nextn2[k2];
      nextn2[k2] = firstn2[ind_1+N*ind_2];
      firstn2[ind_1+N*ind_2] = k2;
      k2 = nextnold;
          if (nprev2 == -1) {
          firstn2[i+N*j] = nextnold;
          }
          else {
          nextn2[nprev2] = nextnold;
          }
      }
      }
    }
    """).build()

div_kernel = prg.div
div_kernel.set_scalar_arg_dtypes([None]*3 + [np.float32, np.uint32] + [None])

L_adj_1_kernel = prg.L_adj_1
L_adj_1_kernel.set_scalar_arg_dtypes([None])

L_adj_2_kernel = prg.L_adj_2
L_adj_2_kernel.set_scalar_arg_dtypes([None])

solve_kernel = solve_prg.solve
solve_kernel.set_scalar_arg_dtypes([None]*7 + [np.float32]*3 + [np.uint32] + [None]*3)

grad_kernel = prg.grad
grad_kernel.set_scalar_arg_dtypes([None]*1 + [np.float32, np.uint32] + [None]*3)

ss_conv_kernel = prg.ss_conv
ss_conv_kernel.set_scalar_arg_dtypes([None]*3 + [np.float32, np.uint32] + [None]*3)

interp_kernel = ib_prg.interp
interp_kernel.set_scalar_arg_dtypes([None]*6 + [np.uint32, np.float32] + [None]*3)
spread_kernel = ib_prg.spread
spread_kernel.set_scalar_arg_dtypes([None]*10 + [np.uint32] + [np.float32]*2 + [np.uint32]*2 + [None]*3)
updaten_kernel = ib_prg.updaten
updaten_kernel.set_scalar_arg_dtypes([None]*6 + [np.uint32] + [np.float32] + [np.uint32]*2)

for ii in range(start_size,end_size+1): 
    #set up grid
    N = 2**ii
    dx = 1.0/N
    rho = 2.0

    x_stag = np.arange(-0.5, 0.5, dx)
    xx_1_stag = np.tile(x_stag, (N, N, 1))
    xx_2_stag = np.tile(x_stag.reshape(N,1), (N, 1, N))
    xx_3_stag = np.tile(x_stag.reshape(N,1,1), (1, N, N))

    x_c = np.arange(-0.5+dx/2, 0.5+dx/2, dx)
    xx_1_c = np.tile(x_c, (N, N, 1))
    xx_2_c = np.tile(x_c.reshape(N,1), (N, 1, N))
    xx_3_c = np.tile(x_c.reshape(N,1,1), (1, N, N))
    
    if PR_INT:
        mesh3 = np.vstack((xx_1_c[np.newaxis],xx_2_c[np.newaxis],xx_3_c[np.newaxis]))

    mu = 0.5
    dt = 0.5*dx
    clock_max = np.int(np.round(.4/dt,0))
    clock_max = 40
    print clock_max

    #set up Lagrangian parameters
    N_r = N/2
    N_s = N/2
    #N_r = 3*N/8
    #N_s = 75*N/16
    b_q = N_r*N_s
    dr = 1./N_r
    ds = 1./N_s
    dtheta = dr*ds
    a_lag = 0.2
    b_lag = 0.25
    g_lag = 0.0625

    #set up Lagrangian grid
    r_lag = np.arange(dr/2, 1, dr)
    s_lag = np.arange(ds/2, 1, ds)
    rr_lag, ss_lag = np.meshgrid(r_lag, s_lag)
    rr_lag = rr_lag.reshape(b_q,1)
    ss_lag = ss_lag.reshape(b_q,1)

    i_lag = np.arange(0, b_q)
    ip_s_lag = np.remainder(i_lag+N_r, b_q)
    im_s_lag = np.remainder(i_lag-N_r, b_q)

    X_1 = (a_lag+g_lag*(rr_lag-0.5))*np.cos(2*pi*ss_lag)
    X_2 = (b_lag+g_lag*(rr_lag-0.5))*np.sin(2*pi*ss_lag)
    X_3 = np.zeros((b_q,1)) 
    #X_3 = np.sin(2*pi*ss_lag) 

    firstn1 = -1*np.ones((N,N))
    nextn1 = -1*np.ones((b_q,1))
    firstn2 = -1*np.ones((N,N))
    nextn2 = -1*np.ones((b_q,1))

    start_time = time.time()
    for k in range(0,b_q):
        s_1_stag = (X_1[k]+0.5)/dx
        s_2_stag = (X_2[k]+0.5)/dx        
        i_1_stag = (np.floor(s_1_stag)).astype(np.uint32)
        i_2_stag = (np.floor(s_2_stag)).astype(np.uint32)
        s_1_c = (X_1[k]+0.5-dx/2)/dx;
        s_2_c = (X_2[k]+0.5-dx/2)/dx;        
        i_1_c = (np.remainder(np.floor(s_1_c),N)).astype(np.uint32);
        i_2_c = (np.remainder(np.floor(s_2_c),N)).astype(np.uint32);
        if (firstn1[i_2_c,i_1_stag]==-1):
            firstn1[i_2_c,i_1_stag]=k 
        else:
            nextn1[k]=firstn1[i_2_c,i_1_stag]
            firstn1[i_2_c,i_1_stag]=k
        if (firstn2[i_2_stag,i_1_c]==-1):
            firstn2[i_2_stag,i_1_c]=k 
        else:
            nextn2[k]=firstn2[i_2_stag,i_1_c]
            firstn2[i_2_stag,i_1_c]=k
    print 'time spent binning points is', time.time()-start_time, 'seconds'

    #move data to gpu

    L_tf, Nyq_filt, p_filt = calc_transf(dx)

    #u_1_soln = 2*np.cos(2*pi*xx_1_stag)*np.sin(4*pi*xx_2_c)
    #u_2_soln = -np.sin(2*pi*xx_1_c)*np.cos(4*pi*xx_2_stag)
    u_1_soln = 2*np.ones((N,N,N))+np.zeros((N,N,N)) #initial conditions
    u_2_soln = np.zeros((N,N,N))
    u_3_soln = np.cos(2*pi*xx_1_c)*np.sin(4*pi*xx_2_c)

    #move data to gpu
    start_time = time.time()
    plan = Plan((N, N, N), queue=queue)
    print 'time spent making plan is', time.time()-start_time, 'seconds'

    u_1_gpu = cl_array.to_device(ctx, queue, u_1_soln.astype(np.float32))
    u_2_gpu = cl_array.to_device(ctx, queue, u_2_soln.astype(np.float32))
    u_3_gpu = cl_array.to_device(ctx, queue, u_3_soln.astype(np.float32))
    test_1_gpu = cl_array.empty_like(u_1_gpu) 
    test_2_gpu = cl_array.empty_like(u_1_gpu) 
    test_3_gpu = cl_array.empty_like(u_1_gpu) 
    u_1_tf_gpu = cl_array.to_device(ctx, queue, u_1_soln.astype(np.complex64))
    u_2_tf_gpu = cl_array.to_device(ctx, queue, u_2_soln.astype(np.complex64))
    u_3_tf_gpu = cl_array.to_device(ctx, queue, u_3_soln.astype(np.complex64))
    L_tf_gpu = cl_array.to_device(ctx, queue, L_tf.astype(np.complex64))
    Nyq_filt_gpu = cl_array.to_device(ctx, queue, Nyq_filt.astype(np.float32))
    p_filt_gpu = cl_array.to_device(ctx, queue, p_filt.astype(np.float32))
    f_1_gpu = cl_array.empty_like(u_1_gpu) 
    f_2_gpu = cl_array.empty_like(u_1_gpu) 
    f_3_gpu = cl_array.empty_like(u_1_gpu) 
    f_1_tf_gpu = cl_array.empty_like(u_1_tf_gpu) 
    f_2_tf_gpu = cl_array.empty_like(u_1_tf_gpu) 
    f_3_tf_gpu = cl_array.empty_like(u_1_tf_gpu) 
    p_gpu = cl_array.empty_like(u_1_gpu) 
    p_tf_gpu = cl_array.empty_like(u_1_tf_gpu) 
    u_1_tf_tilde_gpu = cl_array.empty_like(u_1_tf_gpu) 
    u_2_tf_tilde_gpu = cl_array.empty_like(u_1_tf_gpu) 
    u_3_tf_tilde_gpu = cl_array.empty_like(u_1_tf_gpu) 
    u_1_tilde_gpu = cl_array.empty_like(u_1_gpu) 
    u_2_tilde_gpu = cl_array.empty_like(u_1_gpu) 
    u_3_tilde_gpu = cl_array.empty_like(u_1_gpu) 
    firstn1_gpu = cl_array.to_device(ctx, queue, firstn1.astype(np.int32))
    nextn1_gpu = cl_array.to_device(ctx, queue, nextn1.astype(np.int32))
    firstn2_gpu = cl_array.to_device(ctx, queue, firstn2.astype(np.int32))
    nextn2_gpu = cl_array.to_device(ctx, queue, nextn2.astype(np.int32))
    X_1_gpu = cl_array.to_device(ctx, queue, X_1.astype(np.float32))
    X_2_gpu = cl_array.to_device(ctx, queue, X_2.astype(np.float32))
    X_3_gpu = cl_array.to_device(ctx, queue, X_3.astype(np.float32))
    X_1_tilde_gpu = cl_array.to_device(ctx, queue, X_1.astype(np.float32))
    X_2_tilde_gpu = cl_array.to_device(ctx, queue, X_2.astype(np.float32))
    X_3_tilde_gpu = cl_array.to_device(ctx, queue, X_3.astype(np.float32))
    rr_lag_gpu = cl_array.to_device(ctx, queue, rr_lag.astype(np.float32))
    ip_s_lag_gpu = cl_array.to_device(ctx, queue, ip_s_lag.astype(np.uint32))
    im_s_lag_gpu = cl_array.to_device(ctx, queue, im_s_lag.astype(np.uint32))
    F_1_gpu = cl_array.empty_like(X_1_gpu) 
    F_2_gpu = cl_array.empty_like(X_1_gpu) 
    F_3_gpu = cl_array.empty_like(X_1_gpu) 
    U_1_gpu = cl_array.empty_like(X_1_gpu) 
    U_2_gpu = cl_array.empty_like(X_1_gpu) 
    U_3_gpu = cl_array.empty_like(X_1_gpu) 
    U_1_soln_gpu = cl_array.empty_like(X_1_gpu) 
    U_2_soln_gpu = cl_array.empty_like(X_1_gpu) 
    U_3_soln_gpu = cl_array.empty_like(X_1_gpu) 

    #solve
    queue.finish()
    start_time = time.time()
    for t_step in range(1, clock_max+1):
        solve(t_step)
        if t_step % 1 == 0:
            print t_step
        if PR_INT:
            if t_step % PR_INT == 0:
                #vis_p = p_gpu.get()
                #write_structured_grid("p-%04d.vts"%t_step, mesh3, 
                #point_data=[("p", vis_p)])

                vis_X_1 = X_1_gpu.get()
                vis_X_2 = X_2_gpu.get()
                vis_X_3 = X_3_gpu.get()
                vis_X = np.hstack((vis_X_1,vis_X_2,vis_X_3)).copy()
                data = [
                ("F_1", F_1_gpu.get().reshape(b_q).copy()),
                ("F_2", F_2_gpu.get().reshape(b_q).copy()),
                ]

                file_name = "X-%04d.vtu"%t_step
                compressor = None

                grid = UnstructuredGrid(
                     (b_q, DataArray("points", vis_X, vector_format=VF_LIST_OF_VECTORS)),
                     cells=np.arange(b_q, dtype=np.uint32),
                     cell_types=np.asarray([VTK_VERTEX] * b_q, dtype=np.uint8))

                for name, field in data:
                   grid.add_pointdata(DataArray(name, field,
                   vector_format=VF_LIST_OF_COMPONENTS))

                from os.path import exists
                if exists(file_name):
                   raise RuntimeError("output file '%s' already exists"
                   % file_name)

                outf = open(file_name, "w")
                AppendedDataXMLGenerator(compressor)(grid).write(outf)
                outf.close()

    queue.finish()
    print 'time spent is', time.time()-start_time, 'seconds'

    u_1_soln = u_1_gpu.get()
    u_2_soln = u_2_gpu.get()
    u_3_soln = u_3_gpu.get()
    #print u_1_soln
    #print u_2_soln
    #calculate error
    if ii != start_size: 
        errinf1, errinf2, errinf3, err11, err12, err13, err21, err22, err23 = err(u_1_soln, u_2_soln, u_3_soln, u_1_o, u_2_o, u_3_o)

    #errinf1, errinf2, errinf3, err11, err12, err13, err21, err22, err23 = err_ex(clock_max, u_1_soln, u_2_soln, u_3_soln)
        print 'u Errinf1 is', errinf1, '\n u Errinf2 is', errinf2, '\n u Errinf3 is', errinf3
        print 'u Err11 is', err11, '\n u Err12 is', err12, '\n u Err13 is', err13
        print 'u Err21 is', err21, '\n u Err22 is', err22, '\n u Err23 is', err23

    u_1_o = u_1_soln
    u_2_o = u_2_soln
    u_3_o = u_3_soln
